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ABSTRACT 

Accurate  stress  and  strain  calculations  at  a  notch  usually  require  a  non-linear  finite 
element  analysis  when  local  yielding  has  occurred.  The  strain  energy  density  hypothesis 
is  a  method  to  predict  these  stress  and  strain  values.  This  method  proposes  that  the 
plastic  strain  energy  density  is  equivalent  to  the  strain  energy  density  found  assuming  the 
material  to  be  entirely  elastic.  This  hypothesis  was  evaluated  using  the  finite  element 
method,  which  was  tested  by  comparing  to  exact  solutions  of  elastic  and  elasto-plastic 
problems,  to  calculate  the  stress  and  strain  field  for  two  notched  plates  of  varying  widths 
under  elasto-plastic  loading.  For  both  geometries,  a  plane  stress  and  plane  strain  analysis 
was  performed. 

The  elasto-plastic  strain  energy  density  from  the  finite  element  method  was  found 
to  be  greater  than  that  predicted  by  this  proposal,  which  in  turn  resulted  in  under- 
predicting  the  local  stresses  and  strains.  This  difference  was  greater  for  the  plane  stress 
condition  than  for  the  plane  strain  condition.  Comparisons  were  also  made  with  notch 
stresses  based  the  Neuber  method.  The  two  methods  appear  to  give  an  upper  and  lower 
bound  to  the  actual  stresses  and  strains.  By  combining  the  results  of  the  strain  energy 
density  method  and  the  Neuber  method,  reasonably  accurate  estimates  of  stress  and  strain 
values  can  be  obtained. 


VI 


TABLE  OF  CONTENTS 

I.      INTRODUCTION 1 

A.  FATIGUE  LIFE 1 

B.  STRESS  CONCENTRATION  FACTOR 2 

H.     THEORETICAL  BACKGROUND 5 

A.  NEUBER'S  MODEL 5 

B.  GLINKA  MODEL:  STRAIN  ENERGY  DENSITY  APPROACH 6 

III.STRAIN  ENERGY 9 

A.  STRAIN  ENERGY  PRINCIPLE 9 

B.  DERIVATION  OF  STRAIN  ENERGY  RELATIONSHIP 10 

C.  STRAIN  ENERGY  DENSITY 1 1 

1.  Strain  Energy  Density  for  the  Elastic  Case 12 

2.  Strain  Energy  Density  for  the  Plastic  Case 12 

IV.  VERIFICATION  OF  FINITE  ELEMENT  METHOD 15 

A.  VERIFICATION  OF  ELASTIC  FINITE  ELEMENT  MODEL 17 

1.  Elliptical  Hole  in  an  Infinite  Plate 17 

2.  Conclusion  of  Elastic  Finite  Element  Model 25 

B.  VERIFICATION  OF  ELASTOPLASTIC  FINITE  ELEMENT  MODEL 25 

1.  Uniaxial  Test  of  Stress-Strain  Curve 25 

2.  Infinite  Plate  with  a  Circular  Hole 26 

3.  Comparison  to  Experimental  Data 41 

V.  NUMERICAL  CALCULATION  OF  STRAIN  ENERGY  DENSITY 51 

A.  INTEGRATION  ROUTINE 51 

B.  RESULTS  OF  NUMERICAL  SCHEME  ON  UNIAXIAL  CASE 53 

VI.  CALCULATIONS  AND  RESULTS  OF  STRAIN  ENERGY  DENSITY 57 

A.  NOTCH  GEOMETRY  AND  MATERIAL  SELECTION 57 

B.  FINITE  ELEMENT  MODELING 59 

C.  RESULTS  OF  FINITE  ELEMENT  ANALYSIS 66 

1.  Computational  Procedures 66 

2.  Plane  Stress  Condition 66 

3.  Plane  Strain  Condition 74 

VH.  STRESS  AND  STRAIN  CALCULATIONS 81 

A.  PLANE  STRESS  CONDITION 81 

1.  Finite  Element  Method  Results 81 

2.  Notch  Root  Stress  and  Strain  Calculations  and  Comparisons 84 

B.  PLANE  STRAIN  CONDITION 90 

1.    Finite  Element  Method  Results 90 

vii 


2.    Notch  Root  Stress  and  Strain  Calculations 90 

VIII.  APPLICATION  OF  RESULTS  TO  FATIGUE  CALCULATIONS 101 

A.  STRAIN  LIFE  CALCULATIONS 101 

B.  RESULTS  OF  STRAIN  LIFE  FATIGUE  PREDICTIONS 102 

IX.  CONCLUSIONS 105 

A.  FINITE  ELEMENT  CALCULATIONS 105 

B.  STRAIN  ENERGY  DENSITY  CALCULATIONS 106 

C.  NOTCH  ROOT  STRESS  AND  STRAIN  CALCULATIONS 107 

D.  IMPACT  ON  FATIGUE  LIFE  PREDICTIONS 107 

E.  RECOMMENDATIONS 108 

REFERENCES 109 

INITIAL  DISTRIBUTION  LIST 1 1 1 


Vlll 


LIST  OF  SYMBOLS 

b  fatigue  strength  exponent 

c  fatigue  ductility  exponent 

8,7  Kronecker  delta 

E  Young's  modulus 

Ep  plastic  modulus  (slope  of  plastic  stress-strain  curve) 

Es  secant  modulus  of  plastic  stress-strain  curve 

E\  secant  modulus  of  total  stress-strain  curve 

e  far-field  strain 

e  local  uniaxial  strain 

£,y  total  strain  components 

e'j  elastic  strain  components 

sfj  plastic  strain  components 

eel  equivalent  total  strain 

eee  equivalent  elastic  strain 

ef  equivalent  plastic  strain 

z'f  fatigue  ductility  coefficient 

FEM  finite  element  method 

G  shear  modulus 

K  strength  coefficient 

K,  theoretical  stress  concentration  factor 

K0  stress  concentration  factor 

Kz  strain  concentration  factor 

X  Lame's  constant 

n  strain  hardening  exponent 

Nf  reversals  till  crack  initiation 

a  local  uniaxial  stress 

ix 


Gij  local  stress  component 

Ge  equivalent  stress 

a0  yield  stress 

a'f  fatigue  strength  coefficient 

S  nominal  stress 


Sij         stress  deviator  tensor 


W0  strain  energy  density 

We  strain  energy  density  based  on  elastic  material  properties 

Wp  strain  energy  density  based  on  elasto-plastic  material  properties 

Wa  local  strain  energy  density 

Ws  farfield  strain  energy  density 


ACKNOWLEDGMENT 

The  author  would  like  to  thank  Prof.  Lindsey  for  his  indispensably  guidance  and 
patience  throughout  this  study.  The  author  also  thanks  his  wife,  Renee,  for  her  support, 
encouragement,  and  understanding  during  this  work. 


XI 


I.      INTRODUCTION 
A.         FATIGUE  LIFE 

Fatigue  phenomena  were  first  considered  in  1850  when  it  was  discovered  that 
railroad  car  axles  failed  after  a  certain  amount  of  time  under  cyclic  loads.  The  concept  of 
endurance  limit  was  introduced,  which  states  that  below  a  minimum  load  cyclic  load,  or 
endurance  limit  Se,  some  materials  will  never  fail  [Ref.  1].  To  prevent  fatigue  failure,  one 
only  has  to  design  the  structure  strong  enough  so  that  the  stresses  remain  below  this 
endurance  limit.  However,  certain  structural  applications  require  low  weight,  and  this 
excessive  design  would  be  unfeasible.  Such  is  the  case  with  aircraft.  All  aircraft  are  put 
under  cyclic  loads  every  flight  -  from  repetitive  takeoffs  and  landings,  basic 
maneuvering,  and  flight  gust  loads.  While  aircraft  are  designed  to  withstand  mechanical 
failure  modes  of  large  scale  yielding  and  sudden  monotonic  fracture,  it  is  unrealistic  to 
design  an  aircraft  that  will  never  experience  fatigue. 

It  is  not  surprising  then  that  the  service  lives  of  U.  S.  Naval  Aircraft  are  based  on 
the  fatigue  life  of  the  aircraft  structure.  Due  to  both  a  safety  factor  and  a  cost  factor,  the 
structural  life  of  every  U.  S.  Naval  aircraft  is  thoroughly  tracked.  Engine  or  aircraft 
system  components  can  be  replaced,  but  when  the  main  structural  components  of  an 
aircraft  fail,  there  are  few  options.  The  aircraft  is  either  taken  out  of  service,  or  major  and 
costly  rework  such  as  replacing  the  wings  on  the  A-6E  Intruder  aircraft  is  required  to 
extend  its  life.    Known  weak  spots  are  constantly  inspected  and  reworked,  if  required,  to 
extend  the  life  of  the  aircraft.  These  critical  areas  of  failure,  or  crack  initiation,  usually 
occur  at  stress  concentrations,  such  as  those  occurring  at  notches  or  rivet  holes.  Current 
fatigue  calculations  can  estimate  the  cyclic  life  of  a  component  given  the  range  of  stresses 
and  strains  it  will  undergo.  However,  for  these  calculations  to  be  accurate,  precise  stress 
and  strain  figures  are  required.  This  is  especially  true  when  the  cyclic  fatigue  calculations 
occur  over  tens  of  thousands  of  cycles,  and  any  error  gets  multiplied  many-fold.  Thus, 
accurate  calculations  of  the  stress  at  a  notch  root  is  the  first  step  in  accurate  fatigue  life 
estimation. 


B.    STRESS  CONCENTRATION  FACTOR 

Stress  concentrations  are  stresses  that  are  locally  higher  than  the  far-field  stress  for 
any  loaded  material.  Stresses  increase  and  decrease  due  to  the  geometry  of  a  part,  and 
when  abrupt  geometric  changes  occur,  such  as  a  notch  or  a  hole,  stresses  can  be  two  to 
three  times  greater  than  the  far-field  stresses.  Given  a  defined  geometry  and  applied 
loads,  a  ratio  of  local  (notch  root)  stress  a  to  the  far-field  stress  S  can  be  determined  and 
remains  constant  as  long  as  o  and  S  remain  in  the  linear  range  of  the  stress-strain  curve. 
For  a  notched  specimen,  the  maximum  stress  will  occur  at  the  notch  root,  and  the  ratio  of 
the  maximum  stress  to  the  nominal  stress  is  the  stress  concentration  factor: 

*»-!  (U) 

Throughout  this  thesis,  K,  is  the  theoretical  stress  concentration  factor  based  on  a 
elastic  model.  Even  though  it  is  labeled  theoretical,  its  value  may  be  obtained  from 
analytic  solutions  or  a  finite  element  analysis.  Regardless,  it  is  a  constant  and  holds  for 
any  stresses  in  the  elastic  range. 

The  difficulty  in  calculating  the  stresses  for  a  geometry  such  as  a  notch  or  a  hole 
arise  when  local  yielding  occurs.  Although  an  analytic  solution  or  a  finite  element 
analysis  can  easily  calculate  the  stress  concentration  factor  that  remains  constant  in  the 
elastic  range,  as  soon  as  plastic  yielding  occurs,  this  stress  concentration  factor  decreases 
in  magnitude  as  the  yield  zone  around  the  notch  increases.  This  may  often  be  the  case, 
since  local  yielding  is  allowed  in  the  design  of  the  aircraft,  and  for  stress  concentrations 
of  two  to  three,  it  is  difficult  to  design  a  structurally  efficient  wing  that  will  prevent  local 
yielding.  Once  the  local  stress  departs  the  linear  range  of  the  stress-strain  curve,  it  thus 
becomes  more  difficult  to  calculate  the  local  stresses  and  strains  at  the  notch  root.  A 
finite  element  analysis  is  possible,  but  this  becomes  expensive  and  time  consuming  when 
used  in  high-cycle  fatigue  calculations  that  involve  tens  of  thousands  of  runs. 

In  1961,  H.  Neuber  [Ref.  2]  derived  a  relationship  for  determining  stresses  and 
strains  at  a  notch  that  has  been  loaded  into  the  plastic  range.  Although  Neuber' s 
derivation  involved  a  notch  loaded  in  antiplane  shear,  it  has  been  widely  applied  to 
general  notch  problems.  By  the  1980's,  the  Neuber  method  had  been  adopted  by  virtually 


all  fatigue  analysts,  including  the  Naval  Air  Systems  Command.  However,  over  the  years 
since  Neuber's  original  paper  and  ever  since  it  has  come  into  prominent  use,  many 
investigators  have  proposed  alternative  means. 

One  proposal  by  Glinka  et  al.  [Ref.  3,  4,  5,  6]  is  based  upon  the  concept  that  the 
strain  energy  density  of  the  material  in  the  yielded  zone  is  virtually  the  same  as  the  strain 
energy  density  considering  the  material  to  be  elastic.  This  is  represented  below  in  Figure 
1.1,  where  We  is  the  strain  energy  density  assuming  an  elastic  material  and  equals  the  area 
under  the  linear  curve,  and  Wp  is  the  strain  energy  density  for  an  elasto-plastic  material 
and  equals  the  area  under  the  nonlinear  curve.  This  conjecture  results  in  being  able  to 
calculate  the  stress  concentration  factor  in  the  elasto-plastic  zone  from  the  strain  energy 
density  of  the  elastic  model.  This  proposal  is  based  on  reasoning  that  the  for  local  plastic 
yielding,  there  is  a  relatively  large  volume  of  material  in  the  elastic  region  surrounding 
the  plastic  zone.  Glinka  continued  to  work  on  his  proposal  through  the  1980's,  and 
published  several  papers  applying  his  model  to  plane  stress  and  plane  strain  problems.  In 


W„= }  a  d£ 


>      e 


Figure  1.1.  Representation  of  Strain  Energy  Density  Equivelance  Concept 


1992  W.  N.  Sharpe,  Jr.,  C.  H.  Yang,  and  R.  L.  Trengoning  [Ref.  7]  evaluated  the  Glinka 
relations  with  that  of  Neuber's  for  various  plane  strain  and  plain  stress  configurations 
using  experimental  data  at  the  notch  root.  Their  conclusions  were  mixed,  and  stated  that 


some  cases  were  better  predicted  by  the  Glinka  model,  while  others  were  accurately 
predicted  by  the  Neuber  model. 

Drawing  upon  advanced  finite  element  techniques,  the  proposal  that  the  strain 
energy  density  in  the  plastic  zone  is  equal  to  that  calculated  on  the  basis  of  an  elastic 
solution  has  been  further  tested,  not  only  at  the  notch  root,  but  throughout  the  plastic  zone 
of  the  model.  Elastic  and  Elasto-plastic  problems  with  closed  form  solutions,  along  with 
previously  published  experimental  data  have  verified  the  finite  element  modeling,  which 
was  then  used  to  calculate  strain  energy  density.  Assessments  are  then  made  of  its  impact 
upon  fatigue  life  calculations  of  aircraft  as  compared  to  the  Neuber  approach. 


II.      THEORETICAL  BACKGROUND 

This  thesis  will  concentrate  on  the  Glinka  model  for  calculating  the  stress 
concentration  in  the  elasto-plastic  range.  The  Neuber  model  is  reviewed  below,  and  will 
be  used  for  comparison  of  results  at  the  notch  root  with  that  of  the  Glinka  model. 

A.         NEUBER'S  MODEL 

Neuber  proposed  that  at  the  notch  root  under  plastic  yielding,  the  elastic  stress 
concentration  factor  was  the  geometric  mean  of  the  stress  concentration  factor  and  strain 
concentration  factor,  as  shown  in  Equation  2.1  [Ref.  2]. 


Kt=4KJCz  (2.1) 


where  K„  =  —    and  Kc  =  — 
S  e 


When  the  far-field  stress  S  is  in  the  linear  rage,  this  can  be  rewritten  as: 

Ki    =  KaKe 

S  e 

'      S  S 
or:  (K1S)2  =  Ege  (2.2) 

There  are  two  unknowns  here,  a  and  £.  To  solve  for  them,  the  stress-strain  constitutive 
relationship  is  required.  One  of  the  most  common  stress-strain  relationships  that  can  be 
applied  is  the  Ramberg-Osgood  curve.  For  a  uniaxial  stress  state,  this  is: 


a 
e  =  —  + 
E 


ra:  2.3) 


where  E,  K,  and  n  are  obtained  from  a  curve  fit  to  the  uniaxial  stress-strain  curve. 
Substituting  Equation  2.3  into  Equation  2.2,  the  Neuber  relation  gives  Equation  2.4 
shown  below. 


{K,S)2       a2       fa} 
E  E         [kJ 


(2.4) 


For  the  above  derivation,  the  stress-strain  relationship  using  the  Ramberg-Osgood 
relationship  is  based  on  true  stress,  a  ,  and  true  strain,  e  .  However,  the  computational 
finite  element  analysis  is  based  on  engineering  stresses  and  strains,  a  and  e.  The 
relationship  between  true  and  engineering  stresses  and  strains,  valid  up  to  necking,  is 
given  below  in  Equation  2.5.  It  can  be  seen  that  for  small  strains  such  as  £  <  0.01,  the 
difference  between  true  and  engineering  strains  is  less  than  0.5%,  and  the  difference 
between  true  and  engineering  stresses  is  less  than  or  equal  to  1.0%.  Therefore,  for  small 
strains,  as  is  the  case  in  this  thesis,  the  true  and  engineering  stresses  and  strains  can  be 
interchanged  without  any  significant  error. 

e  =  ln(l  +  e)  a  =  a(l  +  e)  (2.5) 

While  Neuber's  rule  has  been  well  established  as  an  engineering  tool  to  calculate 
notch  stresses  and  strains,  it  has  been  shown  to  overestimate  these  values  [Ref.  3].  The 
accuracy  of  strain  estimation  is  critical  for  fatigue  calculations,  therefore,  a  more  accurate 
method  would  prove  very  beneficial. 

B.         GLINKA  MODEL:  STRAIN  ENERGY  DENSITY  APPROACH 

Glinka  proposed  that  the  energy  density  at  the  notch  root,  calculated  on  the  basis 
of  elasto-plastic  constitutive  laws,  is  equal  to  that  based  on  linear  elastic  constitutive  laws 
for  equivalent  far-field  loading.  In  Chapter  III,  a  detailed  derivation  of  the  strain  energy 
density  will  be  shown,  for  now  it  is  given  below  for  the  uniaxial  case: 


r0=£aJe  (2.6) 


Using  the  above  definition  of  strain  energy  density,  the  strain  energy  at  the  notch  and  far- 
field  regions  are  calculated  (using  a  linear  elastic  stress-strain  relationship  a  =  Ee  ).  For 
the  notch  root,  this  becomes: 


W. 


=  \£Eede  (2.7) 

Jo 


Wa  =  E^-  (2.8) 

2 


w  = 


IE 


(2.9) 


For  the  far-field  region,  g  =  S  and  e  =  e: 


WS=E 


(2.10) 


Ws  = 


IE 


(2.11) 


Rewriting  the  stresses  in  terms  of  strain  energy  density  and  substituting  these  values  into 
the  theoretical  stress  concentration  factor  in  Equation  1.1,  one  gets: 


W 
K,  =  pa 


(2.12) 


However,  Glinka's  hypothesis  is  that  the  strain  energy  density  at  the  root  will  result  in  the 
same  value,  regardless  if  calculated  for  a  linear  elastic  or  an  elasto-plastic  material. 
Therefore,  this  ratio  remains  constant,  even  when  local  yielding  occurs  at  the  notch  root. 
The  argument  is  that  if  the  area  of  local  yielding  is  small,  and  is  surrounded  by  a  large 
volume  of  elastic  material,  then  the  energy  distribution  does  not  change  significantly, 
even  when  local  yielding  occurs. 

For  a  nonlinear  stress-strain  relationship,  Wa  is  found  by  manipulating  the 
integrand  of  Equation  2.6: 

ade  =  e  dG  +  o  dz-z  da  =  d  (a  e)  -  e  da  (2.13) 

Substituting  Equation  2.13  into  Equation  2.6  results  in: 

Wa=j[d(Ge)-edG]  =  Ge-\%dG  (2.14) 

Substituting  the  Ramberg-Osgood  stress-strain  relationship  (Equation  2.3)  into 
Equation  2.14  results  in  the  strain  energy  density  in  terms  of  the  local  uniaxial  stress: 


g2   fn 

w=  —  +  — 

°       E      [k) 


\°f-j 


—  + 
E 


K 


dG 


u;       °2      (X 
°       E      v 


Mi 


1 1  W" " 


_n 


IE     \K 


(ft 


vi  +  l 


(2.15) 


(2.16) 


w„  = 


IE 


\\  +  n 


(a} 


(2.17) 


Substituting  2.17  back  into  the  ratio  for  Kt  (Equation  2.12)  one  obtains  an  equation  in 
terms  of  the  far-field  stress  S,  material  properties  E,  n,  and  K,  and  the  theoretical  stress 
concentration  factor  Kt  that  can  be  solve  numerically  for  the  local  stress  a: 


(KtS)2     g2  ,    2a  (ay 


+ 

E      \  +  n 


(2.18) 


This  expression  applies  only  for  the  uniaxial  case.  In  Chapter  III,  the  strain  energy 
density  is  expressed  in  terms  of  the  general  stress-strain  equations. 

By  inspecting  Equations  2.18  from  Glinka  and  Equation  2.4  from  Neuber,  one  can 
observe  that  the  only  difference  is  the  factor  of  2/(  1+n)  in  the  strain  energy  density 
model.  Since  n  is  less  than  1,  this  term  is  greater  than  1,  and  for  the  left  side  of  these  two 
equations  to  be  equal,  the  local  stress  in  Equation  2.18  must  be  less  than  the  local  stress  in 
Equation  2.4.  Likewise,  if  the  local  stress  in  the  strain  energy  model  are  lower  than  those 
in  the  Neuber  model,  so  will  be  the  local  strains.  In  fact,  Glinka  states  that  the  Neuber 
model  has  been  shown  to  overestimate  the  local  stresses  and  strains  [Ref.  3]. 


III.      STRAIN  ENERGY 

A.         STRAIN  ENERGY  PRINCIPLE 

For  a  general  deformable  body,  external  applied  loads  will  cause  internal  stresses 
and  strains  to  develop  until  an  equilibrium  point  is  reached.  The  result  is  an  internal 
potential  energy,  or  stored  energy  that  is  able  to  do  work.  In  other  words,  potential  strain 
energy  is  the  potential  energy  due  to  internal  stresses  referenced  from  a  zero  stress  state  of 
a  deformable  body.  For  a  system  with  no  losses,  potential  energy  equals  the  work  put 
into  the  system.  To  develop  the  strain  energy  relationship,  consider  a  general  body  with 
externally  applied  loads  as  shown  in  Figure  3.1  below. 

T 


Figure  3.1.  External  Loads  to  a  General  Deformable  Body. 

Given  a  traction  force  T  acting  on  a  portion  of  the  body  dS  with  a  outward  normal 
n  and  that  displaces  5u,  then  the  work  increment  is: 

bW  =  jj(f»Su)dS  (3.1) 

5 

To  develop  potential  strain  energy,  the  basic  governing  equations  are  shown 
below: 

Boundary  Conditions: 

Tx  =cxxnx+axyny+aanz 

Ty=^xynx+<5,yny+(5yznz  (3.2) 


Equations  of  Equilibrium: 


da       da       da 
dx         dy         dz 

3cr       da  vv     dayz 

—^  +  —^  +  —^  =  0  (3.3) 

ox         dy         dz 


dx         dy         dz 
Strain-Displacement  Relationships  (Engineering  Strains): 


du 
£xx=~dx~ 

dv 

Zyy  ~  dy 

dw 

dv     du 
dx     dy 

dw     du 

exz=  —  +  — 

dx     dz 

dw     dv 
dy      dz 

(3.4) 


B.         DERIVATION  OF  STRAIN  ENERGY  RELATIONSHIP 

From  this  point  on,  tensor  notation  will  be  used  as  a  shorthand  notation  to  develop 
the  strain  energy  terms.  The  governing  equations  are  rewritten  below.  The  only  change 
is  in  the  shear  strain  terms,  in  which  the  tensor  shear  strains  are  1/2  of  the  value  of  the 
engineering  shear  strains. 

Boundary  Conditions:  Tt  =  aijnj  ( 3.5 ) 

Equilibrium  Equations:  — —  =  a.  -  =  0  (  3.6  ) 

dx/         J,J 

Strain  Displacement  Relationships:         eiy  =  —  [ul ;j  +  w;,  J  (  3.7  ) 

Returning  to  equation  3.1,  we  substitute  in  the  boundary  conditions  and  expand: 

dW  =  j$Tl§uldS 

s 

dW  =  jj(ai]n])buldS 

s 

dW=\\(aiJdul)njdS 
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At  this  point,  the  Gauss-Divergence  theorem  can  be  applied,  and  the  internal 
strain  energy  can  be  obtained  from  the  externally  applied  loads: 

8W  =  JH(oii8ui)  dV 

v 

V 

When  the  equilibrium  equations  are  substituted  into  the  above  equation,  the  first 
term  goes  to  zero.  The  result  is  expanded  to  form  a  symmetric  and  a  skew-symmetric 
matrix: 

v 
5 w  =  JIf a*  [*  K;  +  K, )  +  2  (&«,;  -  bi„ )]  dV 

V 
V 

where  co  is  defined  as:  co  =  —  \utj  -  uj4  J 

However,  since  co  is  a  skew-symmetric  matrix,  and  since  a  skew-symmetric 
matrix  multiplied  by  a  symmetric  matrix  is  zero,  the  last  terms  disappear  in  the  above 
equations.  The  final  result  is: 

6W  =  JJJayoeiy^  (3-8) 

v 

C.        STRAIN  ENERGY  DENSITY 

Equation  3.8  is  the  differential  strain  energy.  To  find  the  total  strain  energy,  the 
integration  along  each  strain  can  be  accomplished  inside  the  volume  integral.  Strain 
energy  density  is  the  strain  energy  per  unit  volume  and  can  be  represented  as  the 
integrand  of  equation  3.8.  The  8  operator  will  be  changed  to  a  differential,  and  when 
considered  from  a  zero  stress  state,  this  becomes: 

W. 


'.=j/<Vfc5  (3.9) 


where  Ey  is  the  strain  value  at  the  final  stress  state. 
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1.  Strain  Energy  Density  for  the  Elastic  Case 

To  calculate  the  strain  energy  density  for  the  elastic  case,  the  elastic  stress-strain 
relationship,  given  below,  is  substituted  into  equation  3.9: 

a,y=2Ge,;  +  te„5,y  (3.10) 

vE  E  [Ui=j 

with  X  given  as:  X  =  -, r-, r  ;   G  =  —, r  ;   and    5„  =  \ 

5  (l  +  v)(l-2v)  2(1+ v)  ,J      [0;i*j 

The  evaluation  of  the  integral  is  straight  forward: 

W/=Jo£"(2G£y+teu5,,  )</£,,  (3.11) 

w:=GEijelj+^-  (3.12) 

w;  =  Geijeij+-f^biJ  (3.13) 

w;=^(2Ge(;+^£,;5y)  (3.14) 

Replacing  the  term  in  brackets  again  by  the  elastic  stress-strain  relationship,  the  strain 
energy  density  can  be  written  as: 

w;=^aiye,y  (3.15) 

Geometrically,  this  is  the  area  under  the  linear  stress-strain  curve  for  the  case  of 
uniaxial  tension. 

2.  Strain  Energy  Density  for  the  Plastic  Case 

To  evaluate  the  strain  energy  density  for  the  plastic  case,  the  strain  energy  density 
relationship,  equation  3.9,  is  separated  into  an  elastic  term  and  plastic  term: 

W.=w:+Wf=jaM+la9<kl  (3.16) 

The  first  term  is  equation  3.15,  and  the  second  term  can  be  manipulated  in  the  same 
manner  as  Equation  2. 14  to  give: 

W  =-g  ee  +g  ep  -\ep da  (3.17) 

O  ^  'J       'J  <J      ']  J        I]  'J 
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The  stress-strain  relationship  for  elasto-plastic  behavior  for  deformation  theory  based  on 
the  Prandtl-Reuss  equations  is  given  as  Equation  3.18  [Ref.  8]. 


e.y  = 


1  v       * 

2G     ,j       E     kk    v 


3   1 

+ 

2  E 


"       3     ,; 


(3.18) 


Where  Es  is  the  secant  modulus  of  the  effective  stress  versus  plastic  strain  curve: 


E=^ 

5       £p 


(3.19) 


where  the  effective  stress,  Ge,  and  effective  plastic  strain  strain,  e/,  are  defined  as: 
Effective  Stress:  <5e  =  Jfs~s~  (  3.20  ) 

Effective  Plastic  Strain  epe  =  ^fejej  ( 3.21  ) 

where  the  stress  deviator  tensor  is:     s„  =  a„  -  — —b„ 

'J  'J  o        U 

Substituting  the  plastic  strain  component  of  the  relationship  given  in  Equation  3.18  into 
Equation  3.17,  the  strain  energy  density  becomes: 


W  =-a  ee ;+a, 

o      2    'J  'J        'J 


IE 


A 


a  -^5 

»       3 


"J 


3    ( 


IE. 


a  -^5 

»       3     " 


da, 


( 3.22  ) 


If  the  Ramberg-Osgood  uniaxial  stress-strain  curve  is  rewritten  in  terms  of  effective  stress 
and  effective  strain,  it  becomes: 


a      { a 

£  =ee  +ep  =— -  +  \  — 


(  3.23  ) 


Substituting  the  plastic  portion  of  the  effective  strain  of  the  Ramberg-Osgood  relationship 
in  Equation  3.23  into  the  definition  of  the  secant  plastic  modulus,  Equation  3.19,  it  can  be 
rewritten  in  terms  of  the  effective  stress  and  the  Ramberg-Osgood  material  constants. 


J_ 
E. 


K) 


(  1  V     hi 

—     cr  n 

[k]     e 


(  3.24  ) 


Substituting  this  back  into  equation  3.22  results  in  the  strain  energy  density  relationship 
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in  terms  of  stresses  only  (based  on  the  Ramsberg-Osgood  stress-strain  curve): 


W„  =-aI,e'+a,j 

O  n.  1)        IJ  I] 


2K" 


{<>.) 


a  -^5 

u       3      u 


-JAW  k-^5l7Wy(3.25) 
2AT"  V  J 


Before  calculating  the  strain  energy  density  from  the  finite  element  results,  the 
finite  element  program  as  applied  to  notched  geomentries  will  be  verified  in  Chapter  IV. 
Following  verification  of  the  FEM  results,  the  strain  energy  density  Equations  3.9,  3.17 
and  3.25  will  be  numerically  integrated  to  calculate  the  actual  strain  energy  density. 
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IV.     VERIFICATION  OF  FINITE  ELEMENT  METHOD 

To  evaluate  the  strain  energy  density  concept  for  calculating  stress  concentrations 
in  the  plastic  zone,  a  finite  element  model  was  used  to  obtain  the  stress  and  strain  data  as 
an  input  to  the  calculations.  However,  verification  of  the  finite  element  data  was 
paramount  to  ensure  that  a  legitimate  analysis  is  performed.  The  finite  element  program 
used  was  I-DEAS  Master  Series    Release  1.3  and  2.1.  Quadrilateral  plane  stress  and 
plain  strain  elements  were  employed,  and  a  non-linear  stress-strain  analysis  was 
performed. 

The  I-DEAS  Master  Series    offers  several  approaches  to  createing  a  finite 
element  model.  The  general  approach  was  to  create  a  three  dimensional  model  to 
represent  the  physical  specimen.  Next,  an  element  mesh  was  created  by  specifying  the 
type  of  element  to  be  used  and  the  number  of  elements  on  each  side  of  the  surface  being 
meshed.  All  plane  stress  and  plane  strain  analyses  throughout  this  thesis  used  an  eight 
node  quadrilateral  element.  These  elements  are  two-dimensional,  with  nodal  degrees  of 
freedom  consisting  of  translation  in  x  and  y  directions,  and  rotation  about  the  z  axis. 
Only  one  face  of  the  model  was  required  to  be  meshed.  A  mesh  refinement  routine  can 
be  used  to  refine  the  element  shapes  to  reduce  distortion  (i.e.,  skewing  and  stretching), 
and  this  routine  was  used  to  refine  the  finite  element  meshes  for  all  analyses.  After 
generation  of  the  mesh,  boundary  conditions  were  applied  to  represent  external  loads  and 
fixed  displacements.  Note  that  boundary  conditions  can  be  applied  at  lines  of  symmetry; 
thus  reducing  a  block  by  one  half  or  one  quarter  of  the  original  geometry,  resulting  in 
significantly  reduced  computation  time. 

Effective  mesh  generation  is  paramount  in  obtaining  a  correct  solution  for  the 
finite  element  method.  A  coarse  mesh  (few  elements,  large  distance  between  nodes)  will 
not  produce  the  correct  solution  while  an  overly  fine  mesh  (many  elements,  small 
distance  between  nodes)  may  result  in  excessive  computational  time  and  storage. 
Additionally,  a  dense  mesh  near  large  stress  gradients  is  not  only  necessary  for  accurate 
solutions,  but  also  to  provide  the  data  for  effective  post  processing  analysis.  Thus,  for 
each  problem,  a  mesh  study  was  performed  to  determine  the  optimal  mesh  size  and 
distribution  to  use.  This  involved  running  several  finite  element  models  with  varying 
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mesh  sizes  and  distributions,  and  determining  at  what  point  further  mesh  refinements  did 
not  effect  the  solution. 

Two  additional  factors  effect  the  non-linear  portion  of  the  analysis:  the  non-linear 
stress-strain  relationship  and  the  iteration  procedure.  For  a  non-linear  stress-strain 
analysis,  the  I-DEAS  Master  Series    program  uses  input  data  points  to  model  the  stress- 
strain  curve.  The  program  assumes  a  constant  stiffness  between  each  point,  breaking  the 
curve  up  into  individual  linear  sections.  A  maximum  of  twenty  points  can  be  entered. 
To  ensure  a  smooth  curve,  all  twenty  points  were  used  to  define  the  stress-strain  curve, 
with  denser  groupings  at  the  highest  curvature.  Two  basic  iteration  procedures  are 
available:  either  the  Newton-Raphson  method  or  the  Modified  Newton-Raphson  method. 
The  Newton-Raphson  method  was  used,  and  although  this  may  mean  more  computational 
time  per  iteration,  it  generally  iterates  at  a  faster  rate  than  the  Modified  Newton-Raphson 
method,  requiring  fewer  iterations.  For  either  method,  a  convergence  requirement  can  be 
set,  specifying  the  minimum  value  a  function  can  change  from  one  iteration  to  the  next 
before  the  solution  is  considered  to  have  converged. 

Two  basic  types  of  comparisons  can  be  made:  the  first  to  an  analytic  solution,  and 
the  second  to  experimental  data.  The  first  type  provides  a  comparison  to  an  exact 
solution,  and  differences  between  the  finite  element  method  and  an  analytic  solution 
should  be  minimal.  The  advantage  of  this  type  of  comparison  is  that  it  will  show  the 
error  of  the  finite  element  method  compared  to  the  governing  equations  of  solid 
mechanics.  The  comparison  of  the  finite  element  method  to  experimental  data  will 
indicate  how  well  the  finite  element  method  models  actual  mechanical  behavior.  The 
critical  area  in  this  case  is  the  elasto-plastic  relationship,  or  how  well  the  Prandtl-Reuss 
equations  represent  this  relationship.  This  includes  modeling  the  effective  stress  and 
effective  strain  with  a  uniaxial  stress-strain  curve. 

Full  stress  field  measurements  in  solid  mechanics  must  rely  on  surface  techniques 
and  assume  plane  stress  or  plane  strain  conditions.  There  are  also  approximations 
necessary  in  obtaining  the  elasto-plastic  measurements,  such  as  reducing  the  data  via  the 
Prandtl-Reuss  equations.  Another  disadvantage  of  using  experimental  data  is  the  percent 
error  introduced  when  making  measurements.  Despite  these  difficulties  and  experimental 


16 


error,  qualitative  experimental  results  are  indispensable,  and  when  executed  properly,  can 
show  were  the  mathematical  model  differs  from  the  actual  physics.  Numerical 
computational  techniques  will  have  no  validity  unless  they  can  be  shown  to  truly 
represent  the  actual  physical  phenomenon.  Therefore,  despite  the  difficulty  in  obtaining 
experimental  data,  an  attempt  has  been  made  to  compare  the  finite  element  method  to 
previously  published  experimental  results. 

A.         VERIFICATION  OF  ELASTIC  FINITE  ELEMENT  MODEL 

1.  Elliptical  Hole  in  an  Infinite  Plate 

The  first  step  in  the  finite  element  modeling  was  to  ensure  the  program  could 
obtain  accurate  results  for  a  notched  type  of  geometry  under  linear  elastic  loading.  An 
elliptical  hole  in  an  infinite  thin  plate  was  used  for  this  comparison.  Durelli  [Ref.  9] 
gives  the  stress  distribution  around  the  elliptical  boundary  and  compares  this  with 
experimental  data.  Brown  [Ref.  10]  gives  the  complete  stress  field  for  any  uniaxial  load. 
The  stress  distributions  for  a  infinite  plate  with  an  elliptical  hole  are  given  in  elliptical 
coordinates,  which  are  shown  in  Figure  4. 1  and  Figure  4.2,  for  an  ellipse  of  a  =  1  and  b  = 
0.5. 


n  =  i?o» 


TV1 


V  r\  =  60e 

"4  =  0.987 

.0.75-4-4  =  0.783     / 


n  =  180s 


t^  3009 


Figure  4.1.  Elliptical  Coordinate  System. 
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%  =  constant  Oj, 


r|  =  constant 

Figure  4.2.  Normal  Stress  Components  on  an  Element  Referenced  to  Elliptical 

Coordinates. 


The  conversion  from  elliptical  to  Cartesian  coordinates  is  given  by: 

x  =  c  cosh  \  cos  rj 
v  =  c  sinh  ^  sin  T| 

where  c  is  the  focal  distance  of  the  ellipse  given  by: 

c2=a2-b2 

The  stress  around  the  circumference  of  the  elliptical  boundary  is: 

„   ze  sinh  2E„  +  e2^  cos  2ti  - 1 

a   =  S  e^° 22— ! — 

cosh2qo  -cos2r) 


(4.1) 


(4.2) 


were  L,o  is  the  value  of  C,  at  the  boundary  of  the  ellipse. 

The  stress  components,  as  shown  by  Brown  [Ref.  10],  throughout  the  body  are  given  as: 


S 


,2^ 


f  cos  2r|  -  cosh  2\0  +  cosh(2^  -  2l3o ) cos  2r| 
cosh  2\  -  cos  2r] 


sinh  2cj0  sinh  2\  -  sinh  2cj  cos  2r| 
(cosh  2^  -  cos  2rj) 

sinh  2£  -  e{2*-~*)     sinh  2^(cosh  2^0  + 1)  -  sinh  2\  cos  2t| 


(4.3) 


cosh  2\  -  cos  2rj 


(cosh  2^  -  cos  2r|) 


and 
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<5t  =  — 


&. 


cos  2r\  +  cosh  2^0  -  cosh(2^  -  2\0 )  cos  2t] 
cosh  2^  -  cos  2rj 


sinh  2^r)  sinh  2£,  -  sinh  2^  cos  2r| 
(cosh  2^  -  cos  2r|)" 

sinh  2^  -  e^0'2^     sinh  2^(cosh  2^0  + 1)  -  sinh  2^  cos  2r| 


(4.4) 


cosh  2^  -  cos  2rj 


(cosh2^-cos2r|) 


Along  the  x  axis,  a^  =  0,  r|  =  0,  therefore  ax  =  a^  and  ay  =  ar].  On  the  y  axis,  a^  =  0, 
rj  =  7t/2,  therefore  ay  =  c^  and  ax  =  an. 

To  approximate  an  infinite  plate,  a  plate  40  inches  by  40  inches  with  a  central 
elliptical  hole  with  a  major  axis  of  1  inch  and  a  minor  axis  of  0.5  inches  was  used.  A  far- 
field  load  of  5  =  1000  psi  was  applied.  Due  to  symmetry,  only  one  quarter  of  the  block 
was  modeled  (20  inches  by  20  inches).  The  boundary  conditions  required  to  accomplish 
this  symmetry  were  zero  displacement  in  the  y  direction  along  the  x  axis,  and  zero 
displacement  in  the  x  direction  along  the  y  axis.  This  is  shown  below  in  Figure  4.3.  A 


20  in 


0.5  inch 


/tv  A  A  A  A  A 


FEM 
Portion 


IX 


TT 


— * 


k-H 


in 


20  in 


-»    x 


Figure  4.3  Layout  of  model  representing  an  elliptical  hole  in  an  infinite  plate. 
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mesh  of  30  by  30  quadrilateral  plane  stress  elements  was  used.  Figures  4.4  and  4.5  below 
show  the  layout  of  the  finite  element  mesh  used. 


Figure  4.4.  Finite  element  mesh  for  elliptical  hole  in  an  infinite  plate. 


Figure  4.5.  Detail  of  finite  element  mesh  around  elliptical  hole. 

Comparison  of  the  finite  element  results  with  that  of  the  analytical  solution  shows 
very  good  agreement  between  the  two.  Figure  4.6  shows  the  tangential  stress  along  the 
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boundary  of  the  ellipse  plotted  against  the  x-axis  location.  The  finite  element  results 
coincide  with  the  analytic  solution  across  the  complete  boundary. 
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Figure  4.6.   Tangential  Stress  G^  Along  Elliptical  Boundary. 
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Figure  4.7  below  shows  the  percent  error  from  the  analytic  solution.  The  error  is 
generally  less  than  1.0%,  even  at  the  point  of  maximum  stress  (x  =  1).  The  greatest  error 
occurs  near  the  point  of  sign  reversal  (x  =  0.65)  and  is  a  maximum  of  7.45%.  However, 
the  stresses  at  this  point  are  on  the  order  of  5  times  less  than  the  applied  far-field  stresses, 
and  the  actual  solution  crosses  zero.  Therefore,  a  better  representation  of  the  amount  of 
error  is  to  normalize  the  difference  between  the  finite  element  data  and  the  analytic 
solution  by  the  far-field  stress  vice  the  analytic  solution  itself,  as  stated  in  Equation  4.5. 
This  is  shown  in  Figure  4.8. 


Normalized  Stress  Error  = 


Finite  Element  Amalytic  Value 


xl00% 


(4.5) 
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Figure  4. 7  Percent  Error  of  tangential  stress  Gn  along  elliptical  boundary. 
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Figure  4.8  Percent  Normalized  Error  of  tangential  stress  Gv  along  elliptical  boundary. 

A  comparison  was  also  made  of  the  stress  distribution  along  T|  =  0  and  r\  =  nil.  Figure 
4.9  below  shows  the  a^  stress  distribution  along  T|  =  0  for  both  the  finite  element 
calculations  and  the  analytic  solution.  Very  good  correlation  is  shown,  with  the 
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maximum  error  at  1.23%  (see  Figure  4.10).  Figure  4.1 1  below  shows  the  an  stress 
distribution  along  r)  =  tc/2  for  both  the  finite  element  calculations  and  the  analytic 
solution.  Again,  very  good  agreement  of  the  finite  element  method  with  the  analytic 
solution  was  achieved.  The  normalized  error  is  shown  in  Figure  4.12. 
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Figure  4. 9.  Axial  stress  (  On  )  along  rj  =  0. 
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Figure  4.10  Finite  Element  Error  of  Axial  Stress  (  a^  )  Along  rj  =  0. 
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Figure  4.11  Axial  Stress  (  crn  )  Along  r\-n/2. 
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2. 


Conclusion  of  Elastic  Finite  Element  Model 


The  comparison  of  the  finite  element  results  with  the  analytic  solution  of  the 
stress  distribution  around  an  elliptical  hole  shows  that  the  IDEAS  finite  element  program 
can  be  used  to  accurately  analyze  notch  stresses  under  linear  elastic  loading.  Even  though 
the  maximum  error  along  the  hole  edge  was  7.45%,  this  occurred  at  a  point  were  the 
stresses  were  only  3%  of  the  far-field  applied  stress.  The  focus  of  these  finite  element 
analyses  was  to  obtain  accurate  stresses  and  strains  in  the  vicinity  of  the  plastic  zone. 
Therefore,  this  error  at  the  low  stress  regions  becomes  insignificant  when  compared  to  the 
magnitude  of  the  applied  loads  and  of  the  peak  stresses. 


B. 


VERIFICATION  OF  ELASTOPLASTIC  FINITE  ELEMENT  MODEL 


1. 


Uniaxial  Test  of  Stress-Strain  Curve 


The  first  elastoplastic  comparison  of  the  I-DEAS    Software  was  on  a  long  narrow 
block  to  verify  the  input  of  the  stress-strain  curve.  A  finite  element  mesh  of  3  elements 
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by  100  elements  with  the  outer  dimensions  1.0  m  by  0.03  m,  shown  in  Figure  4.13.  This 
simple  case  illustrates  that  the  elasto-plastic  solution  converges  to  the  input  stress-strain 
curve.  Figure  4. 14  shows  the  stress-strain  curve  used,  along  with  the  FEM  calculations. 


Figure  4.13.  Finite  Element  Model  of  Uniaxial  Tension. 
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Figure  4. 14.  Results  of  Uniaxial  Test  Case. 


2. 


Infinite  Plate  with  a  Circular  Hole 


a.         Analytic  Solution  -  Elasto-plastic  Case 


The  elastic  solution  of  the  stress  distribution  of  a  circular  hole  in  an 
infinite  plate  with  a  far-field  applied  load  of  ar  =  a,  as  shown  in  Figure  4. 15,  can  be 
obtained  in  closed  form  and  is  shown  in  Equation  4.6.  To  obtain  the  elasto-plastic  stesses 
and  strains,  both  incremental  and  deformation  theory  have  been  applied  to  this  problem 


( 


.2\ 


G=G 


1- 


ae=a 


<      a2^ 
v      r  J 


(4.6) 
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[Ref.  1 1,  12,  13,  and  14].  Tuba  [Ref.  12]  and  Chakrabarty  [Ref.  14]  show  that  both  these 
solutions  are  identical.  The  solution  below  is  taken  from  Chakrabarty  [Ref.  14]. 


>  a 


Figure  4.15  Infinite  plate  with  circular  hole 


Applying  the  governing  equations  in  polar  form  for  axisymmetric  loading, 
one  obtains  the  stress  components  in  terms  of  constant  ci  and  c?,  that  can  be  solved  for  the 
imposed  boundary  conditions: 


a  =  2c, +^- 


ae=2c2--1 


(4.7) 


The  elastic  solution  shows  that  at  the  hole  edge,  Gq  is  twice  the  far-field  value  a,  with  all 
other  components  being  zero.  From  the  von  Mises  yield  criteria,  Equation  4.8,  yielding 
will  occur  once  G  is  greater  than  one  half  of  the  yield  stress  Go,  and  a  region  of  plastic 
deformation  will  extend  from  the  edge  of  the  hole  out  to  a  certain  radius  of  r  >  a.  If  the 
radius  of  the  elastic/plastic  boundary  is  given  as  c  then  the  boundary  conditions  for 
Equation  4.7  are  Gr  =  Grc  at  r  =  c  and  Gr  =  G  at  r  =  <*>,  where  Grc  is  the  yet  to  be  determined 
value  of  the  radial  stress  at  the  elastic/plastic  boundary. 

G02=G*-GQGr+Gr2  (4.8) 
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The  constant  C3  in  Equation  4.7  can  now  be  solved  for  in  terms  of  Go  and  c,  and  the 
stresses  in  the  elastic  region  reduce  to: 

-2        I  _2 


G=G- 


Si 


1- 


Ga=G  + 


fir* 


1- 


(4.9) 


In  the  plastic  region,  the  Von  Mises  stresses  can  be  expressed  in  parametric  form  as: 


a'  =  s 


aesin<]) 


aa  = 


fi 


CT   COS 


(k       ) 

^6     ; 


(4.10) 


The  stress-strain  constitutive  relationships  given  in  Equation  3.18  are  restated  in  polar 
form: 


1  (  \     3 
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1       1 


kK    e 


£fl  =-(ofl-vOr)  +  - 


^-3(^+^9) 


^e-^^.+^e) 


(4.11) 


The  plastic  secant  modulus  has  been  rewritten  as  follows: 


J_ 


e  -  e 


J_ 
El 


a         a 

j  e  e  j 

were  £/  is  the  secant  modulus  of  the  total  stress-strain  curve.  The  solution  of  the  stresses 
is  independent  of  the  elastic  Poisson's  ratio,  v,  and  therefore  this  can  be  set  equal  to  0.5. 
This  transform  simplifies  the  constitutive  relations,  while  allowing  solution  for  the 
correct  stress  values  [Ref.  13].  With  this  simplification,  and  the  substitution  of  the 
definition  of  the  secant  modulus  as  the  ratio  of  the  effective  stress  to  the  effective  total 
strain,  the  stress-strain  relations  in  the  plastic  region  can  be  expressed  as: 


er  =eesin|(j)-- 


71 


ee  =  e  cos(|> 


(4.12) 


where  er  and  89  are  the  strain  values  for  an  incompressible  material  (v  =  0.5).  Once  the 

stresses  have  been  found,  then  the  actual  strains  can  be  found  using  Equation  4. 1 1  and  the 
correct  value  for  Poisson's  ratio  v. 
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The  strain  ee  is  related  to  the  equivalent  stress  <Je  by  Equation  3.23,  or  a  similar 
stress-strain  curve.  Another  such  empirical  relationship  that  can  be  used  is  a  modified 
Ludwick  curve  [Ref.  14],  which  is  separated  into  an  elastic  and  a  plastic  range: 

£o 
E 

a 


a  =  i 


Ee 


£< 


^ 


vaoy 


(4.13) 


e>^ 


A  plot  of  a  family  of  curves  of  this  form  is  shown  below  in  Figure  4. 16,  where  <Jo  = 
30,000,  E  =  30.0  x  106,  and  n  varies  from  0  to  0.5.     At  the  elastic/plastic  boundary,  the 
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E 

Figure  4.16  Modified  Ludwick  Stress-Strain  Curve. 


auxiliary  angle  $  is  found  by  setting  the  radial  component  of  the  stress  in  Equation  4.9 
equal  with  the  radial  stress  component  in  Equation  4. 10.  The  substitution  of  r  =  c,  and 


a„  =  an  results  in: 


(f)c  =  sin" 


( 


V3a     1 


2a0     2\    a02  j 


(4.14) 


When  stresses  in  Equation  4.10  are  substituted  into  the  equilibrium  condition,  the  result 
is: 
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T^tan(t)  +  v1  =  —  V3-tan(J)  (4.15) 

Ge    dr  dr      2r v  ' 

Likewise,  when  the  strains  in  Equation  4.12,  along  with  the  stress-strain  law  of 
Equation  4.13,  are  substituted  into  the  compatibility  equation,  the  result  is: 

_£-  =  „tan<|)-X  +  _  V3tan<}>-3  (4.16) 

ae    dr  dr      2r v  ' 

The  relationship  between  the  auxiliary  angle  <j) ,  r,  and  n  (material  constant  from 
Equation  4. 13)  can  now  be  found  by  eliminating  the  equivalent  stress  from  Equations 
4.15  and  4.16.  The  resulting  differential  equation  is: 

d§      (V3  -  tan  §)(l  +  V3  n  tan  (})) 
f*  "^r    —  7  ^    \  (4.1/  ) 

dr  2(l  +  ntan"(j)) 

The  solution  to  this  can  be  solved  by  separation  of  variables,  and  can  be  reduced  to  r  as 
function  of  (|>: 


\-n< 


2  1  ~4n  /3  0 

—  =  -r(^/3cos<l)-sin(l))(cos<{)  +  V3ncos<l))1+3',2^     1+3"2  j  (4.18) 

f       a/3 

The  equivalent  stress  must  also  be  found,  and  this  is  obtained  by  eliminating  the  radius 
from  Equations  4.15  and  4.16.  The  resulting  differential  equation  leads  to  the  solution  of 
the  equivalent  stresses  as  a  function  of  <j),  n,  and  a: 

1   dae         (   V3-tan(|) 


=  —n 


(4.19) 


:"('-"). 


ae   3(j>  ^  1  +  V3  n  tan  (j) 

The  relationship  between  (j)  and  Ge  is  now  found  as: 

.  r  rn(l  +  3n)/(l  +  3n2)    \-&TT**\ 

ae  =  /<:ro(cos(})  +  V3nsin(j))  n        e[      1+3n    J  (4.20) 

The  stress  concentration  factor  is  obtained  by  setting  (J>  =  <J)C  and  thus  a(<}))  =  Go  in 
Equation  4.20: 

n(l  +  3n) 


*-?l-«£~W?l  "  "^     <«> 
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For  a  given  applied  load  of  ar  =  a,  the  plastic  stresses  in  Equation  4. 10  can 
be  obtained  by  first  solving  for  tyc  and  K,  as  given  in  Equations  4. 14  and  4.2 1 .  The 
auxiliary  angle  <j)  in  Equation  4.18  is  solved  for  a  given  radius  by  an  iterative  method. 
Finally  <5e  is  solved  from  Equation  4.20,  which  is  then  substituted  into  Equation  4.10  to 
solve  for  the  component  stresses.    For  a  work-hardening  material  (  n  t-  0  ),  the  strains  are 
then  found  by  Equation  4. 1 1. 

b.         Analytic  Calculations 

Two  cases  were  solved  analytically  to  be  used  with  comparison  to  the 
finite  element  solution.  The  first  involved  an  elastic  perfectly  plastic  material  that 
matched  that  used  by  Davis  [Ref.  11]  for  his  numerical  results  using  incremental  theory. 
The  second  case  used  the  same  elastic  properties  but  changed  the  plastic  property  to  a 
work  hardening  material.  Table  4-1  below  shows  the  material  properties: 


Material  Constant 

Value 

Go 

30,000 

E 

30  x  106 

V 

0.3 

a 

22,500 

n,  Case  1 

0 

n,  Case  2 

0.25 

Table  4-1  -  Material  Properties  of  Infinite  Plate  with  hole 

For  the  two  separate  cases,  the  calculated  values  of  <|>c,  c,  and  K,  are  listed  below. 
Figure  4. 17  shows  a  plot  of  Equation  4. 18  for  the  two  cases.  Note  that  there  are  two 
solutions  shown,  and  care  must  be  taken  to  ensure  that  the  correct  root  is  solved  for.  The 
mathematical  software  Maple  V®  was  used  to  solve  (()  for  each  input  radius.  This  program 
allows  the  analyst  to  limit  the  range  in  which  to  search  for  a  solution,  and  therefore  only 
solutions  for  0  <  <j><  60°  were  obtained  without  any  difficulty.  Once  §  was  found  for  a 
specific  radius,  the  calculations  of  Ge,  followed  by  ar  and  Og ,  was  straight  forward. 
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Maple  V    was  programmed  to  make  all  these  calculations  at  a  radial  distance 
corresponding  to  each  node  location  along  the  x-axis. 


Variables 

Case  1  (n  =  0) 

Case  2  (  n  =  0.25  ) 

<j)c 

0.32446 

0.32446 

c 

1.51550 

1.47886 

K, 

1.3333 

1.50195 

Table  4-2  Solution  Variables  for  Case  1  and  Case  2. 


Figure  4.17  Solution  of  Equation  4.18. 


c. 


Finite  Element  Solution 


The  infinite  plate  with  a  circular  hole  was  modeled  using  a  60  degree 
section  on  an  annulus  with  a  inner  hole  radius  of  a  =  1 .0  and  an  outer  radius  of  30.  Due 
to  the  axisymmetric  loading,  a  variety  of  finite  elements  are  available  that  would 
efficiently  model  the  axisymmetric  problem;  however,  the  above  geometry  with 
quadrilateral  plane  stress  elements  was  used  in  order  to  verify  these  elements  and  the 
I-DEAS™  software  for  more  general  notched  geometries. 

The  first  step  was  to  determine  beyond  what  mesh  sizing  the  solution  was 
mesh  independent.  Three  separate  meshes  as  listed  in  Table  4-3  were  analyzed  and  are 
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shown  in  Figures  4. 18  and  4. 19.  The  coarse  and  medium  element  meshes  each  used 
equal  numbers  of  elements  on  opposite  sides  of  the  geometry,  while  the  fine  mesh  had 
one  row  of  elements  at  the  outermost  radius  that  split  into  a  finer  meshing  using 
triangular  elements.  This  was  possible  since  both  the  radial  and  cicumferential  stresses 
rapidly  approach  the  far-field  value,  and  for  the  load  used,  reaches  within  2%  of  the  far- 
field  value  at  a  radial  distance  of  7.5a  and  within  0.5%  of  the  far-field  value  at  15.13a. 
The  finite  element  model  extends  radially  to  30a,  hence  there  is  little  gradient  in  the 
stresses  and  strains  beyond  a  radial  distance  of  15a. 


Mesh  Elements  Nodes 


Coarse 

200 

661 

Normal 

800 

2,521 

Fine 

1378 

4,247 

Table  4-3.  Mesh  sizes  for  hole  in  infinite  plate. 
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Coarse  Mesh 
200  Elements 


Medium  Mesh 
800  elements 


Fine  Mesh 
1378  elements 


Figure  4. 18.  Layout  of  Element  Meshes  for  Hole  in  an  Infinite  Plate. 
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Coarse  Mesh 
200  Elements 


Figure  4.19.  Exploded  Views  of  Element  Meshes  Near  Hole  Boundary. 

The  results  of  the  mesh  comparisons  is  shown  in  Figure  4.20.  The  mesh 
refinement  results  in  only  a  slight  change  in  the  solution,  particularly  at  the  elastic/elasto- 
plastic  boundary.  The  solution  change  from  the  medium  to  fine  mesh  is  insignificant,  and 
it  can  be  concluded  that  the  20  by  40  element  mesh  provides  a  result  that  will  not  improve 
significantly  with  an  increase  in  the  mesh  density.  For  the  rest  of  the  comparisons  to  the 
analytic  solution,  only  the  20  by  40  element  mesh  will  be  used. 
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Figure  4.20  Circumferential  stress  for  various  mesh  sizes. 

The  results  of  the  FEM  were  taken  along  the  one  line  of  elements  at  0  =  0. 
As  would  be  suspected,  the  angular  change  in  the  FEM  solution  was  insignificant,  as 
shown  in  Figure  4.21,  which  shows  the  effective  stress  distribution  throughout  the  finite 
element  model.  The  radius  of  the  elastic/plastic  boundary,  c,  corresponding  to  an 
effective  stress  of  30,000  psi,  was  interpolated  between  nodal  points  to  be  1.5185  inches 
for  case  one,  and  1 .5 107  for  case  two,  with  the  analytic  values  listed  in  Table  4-2.  This  is 
an  error  of  only  0.2%  for  the  former  case,  and  0.59%  for  the  later.  Note  that  the  finite 
element  value  is  linearly  interpolated,  and  though  the  nodal  values  may  be  correct,  this  in 
itself  will  introduce  a  slight  error. 

For  both  case  one  and  case  two,  the  finite  element  solution  was  found  to 
match  the  analytic  solution  very  well.  Figure  4.22  shows  the  finite  element  stresses  with 
the  analytic  curves  superimposed  for  the  perfectly  plastic  model  (case  1).  The  normalized 
error  was  less  than  1%  for  each  data  point  (see  Figure  4.23).  The  stress  distributions  for 
case  2  are  shown  in  Figure  4.24.  Again,  very  good  agreement  is  made  with  the  analytic 
solution,  and  the  normalized  error  remains  less  than  1  %  for  the  work  hardening  material 
(see  Figure  4.25).  For  the  elastic/perfectly-plastic  model,  the  strains  in  the  plastic  region 
cannot  be  obtained  analytically.  However,  for  case  2,  one  simply  solves  Equation  4.11. 
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The  strain  distribution  for  case  2  is  shown  in  Figure  4.26.  As  in  the  case  of  the  stresses, 
the  finite  element  method  produces  very  accurate  results.  The  error  was  normalized  with 
respect  to  {S/E)  as  shown  in  Equation  4.22  for  the  same  reason  the  stress  error  was 
normalized  in  the  preceding  section. 


Normalized  Strain  Error  = 


F  —  £ 

Finite  Element  Analytic  Value 


(S/E) 


xl00% 


(4.22) 


Figure  4.21.  Effective  Stress  for  Case  2. 
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Figure  4.22.   Comparison  of  FEM  stresses  with  analytic  solution  for  case  1. 
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Figure  4.23.  Normalized  error  of  FEM  stresses  versus  analytic  solution  for  case  1. 
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Figure  4.24.   Comparison  ofFEM  stresses  with  analytic  solution  for  case  2. 
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Figure  4.25.  Normalized  error  ofFEM  stresses  versus  analytic  solution  for  case  2. 
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Figure  4.26.   Comparison  of  FEM  strains  with  analytic  solution  for  case  2. 
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Figure  4.27.  Normalized  error  of  FEM  strains  versus  analytic  solution  for  case  2. 
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d.  Conclusion  on  the  Accuracy  of  the  Finite  Element  Method 

Comparisons  of  the  I-DEAS  Master  Series    finite  element  program  to  the 
specific  problems  shown  reveal  the  excellent  results  obtainable  for  both  linear  and  non- 
linear analysis  by  the  finite  element  method.  The  differences  between  the  finite  element 
calculations  with  exact  analytic  solutions  for  both  stress  and  strain  calculations  were 
minimal,  even  in  regions  of  high  gradients.  It  was  also  shown  that  a  mesh-independent 
solution  was  readily  obtainable  after  one  or  two  mesh  refinements.  From  these  results  of 
the  application  of  the  I-DEAS  Master  Series    program  to  model  problems,  follow-on 
finite  element  analyses  for  notch  geometries  were  performed  with  a  very  high  degree  of 
confidence  in  the  accuracy  of  the  solutions. 

3.  Comparison  to  Experimental  Data 

a.  Finite  Element  Modeling  of  Experimental  Specimen 

Although  the  availability  of  recently  published  full-field  stress  data  for 
basic  notched  and  hole  specimens  under  tensile  loads  is  slim  to  non-existent,  several 
works  were  completed  more  than  thirty  years  ago  using  photoelastic  techniques.  Two 
notable  works  are  from  Durelli  and  Sciammarella  [Ref.  15]  and  Theocaris  and  Marketos 
[Ref.  17].  A  comparison  to  the  Theocaris  and  Marketos  experiment  was  chosen  vice  the 
Durelli  and  Sciammarella  due  to  the  fact  that  Theocaris  and  Marketos  were  able  to  show 
the  maximum  peak  tensile  stresses  progressing  away  from  the  notch  root.  This  peak  ay 
stress  progression  occurs  due  to  a  multi-axial  stress  state  and  given  stress  distributions 
allowing  for  a  higher  ay  stress  before  yielding  occurs.  This  was  also  the  result  of  all  the 
finite  element  analyses  completed.  The  test  specimens  used  by  Theocaris  and  Marketos 
were  two  sheets  of  aluminum  alloy  57S,  one  with  a  hole  diameter  to  width  ratio  of  1/2, 
the  other  1/3.  The  later  ratio  of  1/3  was  chosen,  and  the  dimensions  are  shown  below  in 
Figure  4.28. 
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Figure  4.28.  Dimensions  of  Plate  with  a  Central  Hole,  [After  Ref.  17]. 

The  stress-strain  curve  of  the  Aluminum  alloy  57S  as  given  by  Theocaris 
and  Marketos  [Ref.  17]  is  shown  in  Figure  4.29,  with  a  curve  fit  to  the  Ramberg-Osgood 
equation  resulting  in  the  coefficients  listed  in  Table  4-4. 


Coefficient 


Value 


E  (kg/mmz) 

7,000 

K  (kg/mm2) 

33.1 

n 

0.048 

Table  4-4  Ramberg-Osgood  Material  Coefficients  for  Aluminum  57S. 
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Figure  4.29.  Equivalent  Stress-Strain  Curve  for  Aluminum  57S,  [After  Ref.  17]. 

The  first  step,  as  in  the  case  of  the  model  of  the  hole  in  an  infinite  plate, 
was  to  determine  at  what  mesh  density  mesh  independence  was  achieved.  Three  separate 
meshes  were  generated  and  analyzed:  a  20  by  20  mesh,  a  30  by  30  mesh,  and  a  40  by  40 
mesh.  These  are  shown  in  Figure  4.30.  Figure  4.31  shows  the  y-component  of  stress  for 
a  load  well  into  the  plastic  range.  The  change  in  the  FEM  stresses  from  the  30  by  30 
mesh  to  the  40  by  40  mesh  are  insignificant,  hence  the  30  by  30  mesh  will  be  used  for  the 
comparison  to  the  experimental  data. 
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Figure  4.30.   Comparison  of  Finite  Element  Mesh  Sizes  for  a  Hole  in  Finite  Width  Strip. 
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Figure  4.31.  Distribution  of Oyfor  Various  Mesh  Sizes. 

Table  4-5  shows  the  applied  loads  used  by  Theocaris  and  Marketos.  The 
first  load  resulted  in  the  initial  onset  of  plastic  deformation,  while  the  last  load  creates  a 
plastic  zone  that  extends  a  distance  3/4  of  the  hole  radius  from  the  edge  of  the  hole  and 
covering  approximately  1/3  of  the  minimum  cross  section. 


Load  Set 

Applied  Force 

(N) 

S  (MPa) 

SNET  (MPa) 

I 

19162.97 

67.89 

101.84 

n 

21384.95 

75.76 

113.65 

m 

25520.79 

90.42 

135.63 

rv 

30399.66 

107.70 

161.55 

V 

34286.34 

121.47 

182.21 

VI 

38926.48 

137.91 

206.87 

Table  4-5.  Load  history  for  Theocaris  Photoelastic  Experiment,  [Ref.  17] . 
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b.  Finite  Element  Results 

A  comparison  of  the  ay  stress  was  made  for  all  six  load  sets,  and  these  are 
shown  in  Figures  4.32  and  4.33.  Several  differences  are  noticeable  between  the  FEM 
results  and  the  Theocaris  and  Marketos  data.  First,  although  Theocaris  and  Marketos 
show  a  peak  ay  stress  that  progresses  inward  as  the  plastic  deformation  increases,  they 
also  show  that  it  initially  decreases  before  reaching  a  maximum  for  the  last  three  applied 
loads.  The  FEM  shows  that  the  maximum  value  of  ay  moves  away  from  the  hole  edge, 
but  it  also  shows  that  ay  constantly  increases  until  it  peaks.  The  second  difference 
between  the  FEM  and  the  Theocaris  and  Marketos  data  is  the  magnitude  of  the  decrease 
in  ay  near  the  edge  opposite  the  hole.  This  results  in  significant  disagreement  between 
the  two  values,  approaching  an  80%  difference  for  the  first  applied  load.  Figure  4.34 
shows  the  difference  between  the  FEM  calculations  and  the  experimental  data.  One  test 
of  the  accuracy  of  both  results  is  to  determine  if  equilibrium  has  been  satisfied  at  the 
minimum  cross  section.  The  stress  distribution  curves  were  numerically  integrated  to 
determine  the  resulting  force.  These  results  are  listed  in  Table  4-6.  The  values  are  only 
for  half  of  the  plate,  hence  the  total  applied  force  will  be  twice  these  values.  Even  though 
these  calculations  are  only  approximate,  it  is  easily  seen  that  the  FEM  has  satisfied 
equilibrium,  while  the  Theocaris  and  Marketos  results  have  underestimated  the  stress 
distribution  for  the  first  three  load  sets.  It  should  also  be  noted  that  the  photoelastic 
analysis  shows  fringes  in  the  regions  of  strain  gradients,  and  since  this  is  a  region  of 
relatively  uniform  stress  and  strain,  the  accuracy  of  the  method  degrades. 
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Figure  4.32.  Distribution  of  Oyfor  Theocaris  model,  load  sets  I  to  III. 
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Figure  4.33.  Distribution  of  Gyfor  Theocaris  model,  load  sets  IV  to  VI. 
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Figure  4.34.  Difference  of  FEM  and  Theocaris  data  for  (7V. 
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Load 

Applied 

FEM  Calculations 

Theocaris  Data 

Set 

Force 

Fy  at  Far 

Fy  at  minimum 

Error 

Fy  at  minimum 

Error 

Edge (N) 

section  (N) 

section  (N) 

I 

19162.97 

19166.13 

0.02% 

17254.83 

-9.96% 

n 

21384.95 

21387.70 

0.01% 

20004.95 

-6.45% 

m 

25520.79 

25522.31 

0.01% 

24975.24 

-2.14% 

IV 

30399.66 

30401.50 

0.01% 

30473.38 

0.24% 

V 

34286.34 

34286.69 

0.00% 

35225.31 

2.74% 

VI 

38926.48 

38925.78 

0.00% 

39660.45 

1.89% 

Table  4-6.  Equilibrium  calculations  at  minimum  cross-section. 

c.  Conclusion  on  the  Accuracy  of  the  Finite  Element  Method 

The  finite  element  analysis  compared  favorably  with  the  experimental  data 
of  Theocaris  and  Marketos.  Within  20  mm  of  the  hole  edge,  the  error  was  less  than  10% 
for  all  but  one  load  set.  Additionally,  the  results  at  the  edge  itself  were  within  4%.  In 
addition  to  the  experimental  errors  referred  to  in  the  beginning  of  this  Chapter,  the  data 
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used  in  this  study  was  taken  from  published  graphs,  resulting  in  an  additional  error  added 
to  the  comparison.  Considering  these  shortcomings,  the  FEM  analysis  provided 
quantitative  results  as  reasonably  as  could  be  expected,  and  matched  the  qualitative  trends 
exceptionally  well. 
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V.     NUMERICAL  CALCULATION  OF  STRAIN  ENERGY  DENSITY 

A.         INTEGRATION  ROUTINE 

In  order  to  analyze  the  strain  energy  density  from  the  finite  element  data,  an 
integration  routine  must  be  used  to  calculate  the  plastic  strain  energy  density  as  shown  in 
equation  3.25.  The  accuracy  of  any  integration  routine  will  be  dependent  on  the  number 
of  data  points  integrated  and  the  scheme  implemented.  This  is  especially  true  at  the 
'knee'  of  the  stress-strain  curve,  or  just  past  yielding  where  the  curve  bends.  To  increase 
accuracy,  a  third  order  accurate  integration  scheme  was  used  by  approximating  the  first 
and  second  derivative  of  each  point  by  a  central  difference  scheme,  and  using  these  values 
to  perform  the  numerical  integration.  This  scheme  is  similar  to  Simpson's  rule,  but 
applied  for  non-equally  spaced  points.  To  develop  the  integration  routine,  a  function  I(x) 
is  assumed  in  which: 

I(x)  =  \f{x)dx  (5.1) 

Expanding  /(x)  about  the  i   point  by  Taylor  series,  we  get: 

Ax  Ax 

/(x,+Ax,)  =  /(x1)+Ax,/'(x1)+— j-r(x,)+— j- /'"(*,  )+tf.ar. 

2  3  (5-2) 

Ax  Ax 

/(x,  -  Ax,_, )  =  /(x, ) -  Ax^IXxh )  +  —£- /"(*, ) Jj±- /'"(x, )  +  H.O. T. 

were  the  subscripts  refer  to  each  data  point  as  shown  in  Figure  5.1. 
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Figure  5. 1  Integration  points  for  numerical  scheme. 


However,  from  Equation  5.1,  we  know  that: 

/'  =  /(*)    r\x)  =  f\X)  rxx)  =  f\x) 

Substituting  these  into  Equation  5.2  and  solving  for  the  integration  between  the  i    -  1  and 
the  i    +  1  point  and  defining  this  value  as  k  results  in: 

/,.=  ]f(x)dx 


I.  =  I(x,  +Axt)-  I(xt  -  Axw )  +  H. O. T. 


/,.  =  (A*,. +  Azw  )/(*,)  + 


2! 


/'(*,)■ 


Ax,3  +  Ax._,3 
3! 


f"(xt)+H.O.T.       (5.3) 


The  first  and  second  derivatives  of  the  function  are  approximated  by  a  central  difference 
scheme  as  follows: 

A*w2(/(*w )  -  /(*, ))  +  Axl\f(xl )  -  /(*,_, )) 


/'(*,)  = 


Ax1_1Ax,.(Ax._1+AxI.) 


(5.4) 


and 


/,(x)_2^l-,(/(^1)-/(xi))  +  Ax,(/(x,_,)-/(xl)) 
Ax^Ax^Ax^  +  Ax,) 


(5.5) 
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Equations  5.4  and  5.5  are  substituted  into  Equation  5.3,  then  summed  over  every  other 
data  point  to  obtain  the  final  integration  value: 


j7(*)=£/,  (5.6) 


.=2 
even  i 


B.         RESULTS  OF  NUMERICAL  SCHEME  ON  UNIAXIAL  CASE 

A  comparisons  of  the  above  numerical  integration  routine  to  a  standard 
trapezoidal  integration  routine  were  made,  using  the  uniaxial  strain  energy  density 
equation  with  the  Ramberg-Osgood  stress-strain  relationship.  The  stress-strain  curve 
integrated  is  shown  below  in  Figure  5.2,  with  the  limits  for  both  stress  and  strain  labeled. 
The  material  constants  were  n  =  0.053,  K  =  95ksi,  and  E  =  10  xlO3  ksi.  For  the  modified 
Simpson's  rule  integration  routine,  two  forms  of  the  strain  energy  integral  were  used, 
Equation  2.6  and  2.14.  The  actual  value  is  given  by  Equation  2.17.  Note  that  for 
integration  by  the  trapezoid  rule,  both  of  these  forms  of  the  strain  energy  density  produce 
identical  values;  however,  for  the  applied  integration  routine,  an  upper  and  lower  bound 
is  calculated.  This  is  shown  in  Figure  5.3,  which  applied  the  routine  to  equally  spaced 
points.  Note  that  the  actual  implementation  of  this  routine  will  produce  improved  results 
due  to  a  more  efficient  distribution  of  the  data  points;  i.e.,  few  points  along  the  linear 
portion  and  a  closer  distribution  in  the  region  of  higher  curvature.  Since  the  modified 
Simpson's  integration  routine  produced  better  results  than  the  trapezoidal  rule  for  higher 
data  points,  this  was  chosen  as  the  means  of  integration  to  calculate  the  strain  energy 
density.  Additionally,  the  basic  form  of  the  strain  energy  density  integral  as  given  in 
Equation  2.5  was  able  to  be  incorporated  within  the  I-DEAS  Master  Series    post- 
processing module;  hence,  this  was  the  final  form  used  to  calculate  the  strain  energy 
density. 
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Figure  5.2.  Stress-Strain  Curve  for  Integration  Comparison. 


820  i 
810 

800  ■] 


I   790 


780 


c   770  H 

3   760 

q 

cd 

a   750 

u 

3   740  <f 


730 

720  ti- 


/ 


\ 


\ 


/ 


\ 


\ 


—  £r- 

Modified  Simpson's  Rule  for  Equation  2.5 

o 

Modified  Simpson's  Rule  for  Equation  2.13 

-□  - 

Trapezoidal  Rule 

Actual  Value 

"*. 


A-  — -  __ 


O 


—      B- 


^rrfr 


A-    —    —   A 


9  11  13  15  17 

Number  of  Equally  Spaced  Stress  Data  Points 


19 


21 


Figure  5.3.   Comparison  of  Integration  Routines. 
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VI.      CALCULATIONS  AND  RESULTS  OF  STRAIN  ENERGY  DENSITY 

In  order  to  test  the  validity  of  the  Glinka  strain  energy  density  proposal,  the  strain 
energy  density  for  a  given  loading  was  calculated  based  on  an  elasto-plastic  material  (Wp) 
and  an  elastic  material  (We),  and  comparisons  between  the  two  calculations  were  made. 
The  strain  energy  density  calculations  were  performed  at  each  node  throughout  the  FEM 
model,  with  the  plastic  strain  energy  density  (Wp)  being  numerically  integrated  at  each 
load  step  as  shown  in  Chapter  V.  In  addition  to  the  strain  energy  density  calculations, 
comparisons  were  made  between  the  notch  stresses  and  strains  based  on  the  finite 
element  method  and  values  obtained  via  both  the  Glinka  and  the  Neuber  method.  These 
comparisons  will  be  made  in  Chapter  VII.  The  finite  element  analysis  was  performed  on 
a  plate  with  symmetrical,  semi-circular  edge  notches. 

A.         NOTCH  GEOMETRY  AND  MATERIAL  SELECTION 

Two  separate  geometries  were  evaluated,  one  had  a  notch  radius  of  1 .0  inch,  and  a 
plate  width  of  6.0  inches  (r/D  =  1/6),  and  the  second  had  a  notch  radius  of  1.0  inch  and  a 
plate  width  of  10.0  inches  (r/D  =  1/10).  These  two  layouts  are  shown  below  in  Figure 

6. 1,  with  the  shaded  regions  representing  the  finite  element  geometry.  For  each 
geometry,  plane  stress  (thin  plate)  and  plane  strain  (thick  plate)  analyses  were  performed. 
The  strain-hardening  material  used  modeled  7075-T6  Aluminum.  This  was  represented 
in  the  I-DEAS  Master  Series    finite  element  program  as  20  data  points,  shown  in  Figure 

6.2.  This  figure  also  shows  the  yield  stress,  determined  to  be  66  ksi  based  on  a  plastic 
strain  offset  of  0.002.  However,  it  should  be  noted  that  the  stress-strain  curve  departs 
from  linearity  at  40  ksi.  The  20  data  points  were  obtained  by  curve  fitting  the  Ramberg- 
Osgood  equation  to  actual  stress-strain  data  from  the  Military  Handbook  V  [Ref.  1 10]; 
then  calculating  these  points  from  the  resulting  Ramberg-Osgood  equation.  Although  this 
may  have  resulted  in  the  stress-strain  curve  used  by  the  FEM  analysis  differing  from  the 
actual  7075 -T6  data,  it  ensured  that  the  notch  root  strains  and  plastic  strain  energy  density 
could  be  calculated  from  the  finite  element  stresses  using  the  Ramberg-Osgood  equation. 
This  enabled  correlating  predicted  strains  based  on  the  Ramberg-Osgood  equation  and 
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either  the  Neuber  or  Glinka  method  to  the  finite  element  results.  The  material  property 
coefficients,  including  values  resulting  from  the  data  fit  to  the  Ramberg-Osgood  equation, 
are  listed  in  Table  6- 1 . 
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Figure  6.1.  Notch  Geometry  and  Finite  Element  Boundary  Conditions. 
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Figure  6.2  Stress-Strain  Curve  for  707 5 -T6  Aluminum. 


Coefficient 

Value 

E  (ksi) 

10,000.0 

K  (ksi) 

92.0 

a0  (ksi) 

66.18 

n 

0.053 

Table  6-1.   7075-T6  Material  Constants. 


B. 


FINITE  ELEMENT  MODELING 


Due  to  the  two  lines  of  symmetry  for  this  geometry,  the  finite  element  model  was 
reduced  to  one  quarter  of  the  physical  model  by  applying  appropriate  constraints  along 
each  boundary,  as  shown  in  Figure  6.1.  For  each  configuration,  21  increments  were  used 
to  increase  the  loading  from  an  initial  nominal  stress  of  12  ksi  to  a  final  value  of  49.5  ksi, 
or  75  percent  of  the  yield  stress.  This  load  increment  ensured  convergence  of  the  finite 
element  solution  and  provided  a  small  enough  step  to  numerically  integrate  the  strain 
energy  density  with  reasonable  accuracy.  The  first  two  load  points  were  in  the  elastic 
range,  and  plastic  deformation  began  at  the  third  load  step.  The  stress  concentration  factor 
K,  was  determined  from  the  finite  element  analysis  at  the  first  load  step,  and  found  to  be 
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2.042  and  2.421  for  the  narrow  and  the  wide  plate,  respectively.  Table  6-2  shows  the 
load  schedule  for  each  of  the  geometries. 

As  in  the  previous  finite  element  analyses,  a  mesh  comparison  was  performed  for 
each  of  the  geometries.  The  number  of  elements  and  nodes  for  each  geometry  used  in  the 
comparison  are  shown  below  in  Table  6-3.  For  the  wide  plate  (r.D  =  1:10),  the  medium 
mesh  was  obtained  directly  from  the  coarse  mesh  by  an  I-DEAS  Master  Series    routine 
of  splitting  elements  in  the  region  of  highest  elastic  strain  energy  density.  The  mesh 
layouts  are  shown  in  Figures  6.3  and  6.4.  Comparisons  of  both  stresses  and  strains  were 
made  for  each  of  the  mesh  sizes  at  load  step  2 1 .  The  comparison  showed  that  the  mesh 
refinement  resulted  in  little  change  in  the  stresses  and  strains  throughout  the  model. 
Figures  6.5  and  6.6  show  the  stress  and  strain  distributions  along  the  minimum  cross 
section  for  the  each  plate.  For  the  narrow  plate  at  the  notch  root  itself,  the  stress 
decreased  by  0. 16%  going  from  the  coarse  to  medium  mesh,  and  decreased  by  0.03% 
going  from  the  medium  to  fine  mesh.  For  the  strains,  the  medium  mesh  resulted  in  a 
0.24%  decrease  in  strain  compared  to  the  coarse  mesh,  and  there  was  no  change 
compared  to  the  fine  mesh.  Similar  results  are  shown  for  the  mesh  refinement  for  the 
wide  plate.  The  mesh  adaptation  used  on  the  coarse  mesh  to  obtain  the  medium  mesh 
resulted  in  a  decrease  of  the  notch  root  stress  by  only  0.22%.  The  strain  at  the  notch  root 
decreased  0.58%  from  the  medium  mesh  to  the  coarse  mesh.  These  results  showed  that 
the  medium  mesh  for  both  geometries  was  sufficient  to  provide  mesh  independent 
solutions.  Hence,  for  all  further  finite  element  analysis,  the  medium  meshes  were  used. 
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Load 
Step 

Notch  r.\D  = 

1:6 

Notch  r:D  =  1 

:10 

P  (ksi) 

S  (ksi) 

Kf  (ksi) 

P  (ksi) 

S  (ksi) 

K,S  (ksi) 

1 

8.00 

12.00 

24.50 

9.60 

12.00 

29.05 

2 

13.00 

19.50 

39.81 

15.60 

19.50 

47.21 

3 

15.00 

22.50 

45.94 

18.00 

22.50 

54.47 

4 

16.00 

24.00 

49.00 

19.20 

24.00 

58.10 

5 

17.00 

25.50 

52.06 

20.40 

25.50 

61.73 

6 

18.00 

27.00 

55.13 

21.60 

27.00 

65.36 

7 

19.00 

28.50 

58.19 

22.80 

28.50 

68.99 

8 

20.00 

30.00 

61.25 

24.00 

30.00 

72.63 

9 

21.00 

31.50 

64.31 

25.20 

31.50 

76.26 

10 

22.00 

33.00 

67.38 

26.40 

33.00 

79.89 

11 

23.00 

34.50 

70.44 

27.60 

34.50 

83.52 

12 

24.00 

36.00 

73.50 

28.80 

36.00 

87.15 

13 

25.00 

37.50 

76.56 

30.00 

37.50 

90.78 

14 

26.00 

39.00 

79.63 

31.20 

39.00 

94.41 

15 

27.00 

40.50 

82.69 

32.40 

40.50 

98.04 

16 

28.00 

42.00 

85.75 

33.60 

42.00 

101.68 

17 

29.00 

43.50 

88.81 

34.80 

43.50 

105.31 

18 

30.00 

45.00 

91.88 

36.00 

45.00 

108.94 

19 

31.00 

46.50 

94.94 

37.20 

46.50 

112.57 

20 

32.00 

48.00 

98.00 

38.40 

48.00 

116.20 

21 

33.00 

49.50 

101.06 

39.60 

49.50 

119.83 

Table  6-2.  Load  Increment  Schedule. 


Geometry 

Mesh 

Element 

Nodes 

Narrow  Plate 

Coarse  Mesh 

298 

965 

(r.D=  1:6) 

Medium  Mesh 

1,128 

3,507 

Fine  Mesh 

2,291 

7,046 

Wide  Plate 

Coarse  Mesh 

456 

1,453 

(r.D=  1:10) 

Medium  Mesh 

1,064 

3,233 

Table  6-3.  Number  of  Elements  and  Nodes  for  Mesh  Comparison. 
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Figure  6.4.  Mesh  Layouts  for  Wide  Plate. 
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Figure  6.5.  Distribution  of  Gyfor  Mesh  Comparison  of  Narrow  Plate  FEM  Model. 
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Figure  6.6.  Distribution  of £yfor  Mesh  Comparison  of  Narrow  Plate  FEM  Model. 
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Figure  6.7.  Distribution  of  Gy  for  Mesh  Comparison  of  Wide  Plate  FEM  Model. 
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Figure  6.8.  Distribution  of  ey  for  Mesh  Comparison  of  Wide  Plate  FEM  Model. 
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C.         RESULTS  OF  FINITE  ELEMENT  ANALYSIS 

1.  Computational  Procedures 

The  elastic  strain  energy  density  (We)  at  each  node  was  initially  calculated  at  the 
first  load  step.  To  determine  We  at  later  load  steps,  this  value  was  simply  multiplied  by  a 
factor  of  (Pj/Pq)  ,  where  Pt  is  the  load  at  the  z'th  step,  and  P0  is  the  initial  load.  Plastic 
strain  energy  density  (Wp)  was  calculated  incrementally  with  the  integration  procedure 
shown  in  Chapter  V.  To  calculate  a  strain  energy  density  increment  AWP,  two  new  data 
points  were  required,  in  addition  to  last  point  of  the  previous  increment.  Since  the  first 
load  step  was  under  elastic  conditions,  Wp  is  equal  to  We;  hence,  the  first  load  step 
resulted  in  the  first  plastic  strain  energy  density  calculation  without  requiring  any 
integration.  Subsequent  calculations  were  performed  at  all  odd  load  steps. 

2.  Plane  Stress  Condition 

For  both  plate  widths,  the  plastic  strain  energy  density  was  found  to  be  greater 
than  the  elastic  strain  energy  density  in  the  vicinity  of  the  notch  root.  Figures  6.9 
and  6. 10  are  plots  of  Wp  at  load  step  21  for  the  narrow  and  wide  plates,  respectively, 
under  plane  stress  conditions.  These  show  that  the  plastic  strain  energy  density 
throughout  the  model  has  its  maximum  value  at  the  notch  root,  but  rapidly  approaches  the 
far-field  value  away  from  the  notch.  Figures  6.1 1  and  6.12  show  Wp  for  several  load  steps 
along  the  minimum  cross  section  of  the  plates,  with  We  included  at  the  final  load  step. 
This  shows  that  not  only  does  We  give  under  estimated  values  at  the  notch  root,  but  also 

W-We 
WEKKOR=-J^-L  (6.1) 

follows  a  different  distribution  shape  than  Wp.  At  the  notch  root  itself  for  the  final  load 
step  (S  =  0.75a0),  the  difference  between  Wp  and  We  reaches  a  value  of  16.8%  and  23.2% 
for  the  narrow  plate  and  wide  plate  respectively.  Figures  6.13  and  6.14  show  contour 
plots  of  the  difference  in  Wp  as  compared  to  We,  as  shown  in  Equation  6. 1 ,  so  that  positive 
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values  indicate  that  the  actual  Wp  is  higher  than  that  estimated  by  the  Glinka  Method. 
These  values  are  taken  at  the  final  load  step  (S  =  0.75  (Jo).  In  the  plastic  region,  the 
greatest  difference  between  the  two  calculations  occurs  not  at  the  notch  root  itself,  as  may 
be  expected,  but  slightly  offset  along  the  x-axis.  This  maximum  error  is  also  on  the  order 
of  twice  the  value  of  that  at  the  notch  root.  The  location  of  the  maximum  difference 
between  Wp  and  We  corresponds  fairly  well  with  the  location  of  the  maximum  ay  value.  It 
should  be  noted  that  although  the  contour  plots  show  the  maximum  global  error  occurring 
slightly  above  the  notch  near  the  plate  edge,  this  is  also  a  region  of  very  low  to  zero  strain 
energy  density,  and  therefore  these  errors  are  actually  insignificant.  Figures  6.15  and  6.16 
show  the  error  in  strain  energy  density  for  each  load  step  across  the  x-axis  (at  y  =  0  ). 
These  plots  show  that,  as  previously  stated,  the  maximum  error  occurs  at  the  notch  root 
for  small  plastic  yielding,  then  gradually  progresses  inward  along  the  x-axis.  Global 
errors  for  the  previous  load  steps  followed  the  save  trends  as  that  of  the  final  load  step. 
Starting  with  zero  error  at  the  initiation  of  plastic  deformation,  the  regions  of  significant 
error  (greater  than  1%)  start  at  the  notch  root,  and  move  inward,  while  at  the  same  time 
proceeding  at  an  angle  toward  the  vertical  centerline  of  the  plate.  This  trend  also 
corresponds  to  the  region  of  high  ae  and  the  growth  of  the  plastic  zone  (this  growth  is 
well  documented  by  Theocaris  and  Marketos  [Ref.  18]). 
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Figure  6.9.  Plastic  Strain  Energy  Density  Wpfor  Narrow  Plate  in  Plane  Stress 

at  Load  Set  21. 
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Figure  6.10.  Plastic  Strain  Energy  Density  Wpfor  Wide  Plate  in  Plane  Stress  at 

Load  Set  21. 
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Figure  6.11.  Strain  Energy  Density  along  x-axisfor  Narrow  Plate  in  Plane  Stress. 
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Figure  6.12.  Strain  Energy  Density  along  x-axisfor  Wide  Plate  in  Plane  Stress. 
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Figure  6.13.  Difference  Between  Wp  and  Wefor  Narrow  Plate  in  Plane  Stress  at 

Load  Step  21. 
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Figure  6. 14.  Difference  Between  Wp  and  Wefor  Wide  Plate  in  Plane  Stress  at 

Load  Step  21. 
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Figure  6. 15.  Percent  Difference  Between  Wp  and  Wefor  Narrow  Plate  in  Plane 
Stress  Along  Minimum  Cross  Section. 
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Figure  6. 16.  Percent  Difference  Between  Wp  and  Wefor  Wide  Plate  in  Plane 
Stress  Along  Minimum  Cross  Section. 
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3.  Plane  Strain  Condition 

The  results  of  the  strain  energy  density  calculations  showed  better  agreement 
between  Wp  and  We  than  for  the  plane  stress  condition.  Figures  6.17  and  6.18  are  plots  of 
Wp  at  the  final  load  step.  As  in  the  plane  stress  case,  there  appears  a  high  peak  value  at 
the  notch  root,  with  a  rapid  drop-off  to  the  far-field  value.  Figures  6. 19  and  6.20  show  Wp 
for  several  load  steps  along  the  minimum  cross  section  of  the  plates,  with  We  included  at 
the  final  load  step.  Closer  agreement  between  the  two  strain  energy  densities  is  shown 
than  for  the  plane  stress  condition.  The  plane  strain  problem  results  in  a  configuration 
that  is  physically  more  constrained  than  that  of  the  plane  stress  problem,  and  hence  the 
amount  of  plastic  growth  at  the  notch  under  plane  strain  conditions  will  be  less  than  that 
of  the  plane  stress  condition.  Therefore,  the  basis  that  the  strain  energy  density 
distribution  in  the  plastic  region  remains  relatively  unchanged  due  to  a  high  volume  of 
elastic  material  surrounding  the  plastic  region  should  be  even  more  valid  for  the  plane 
strain  condition  than  that  of  the  plane  stress  condition.  This  was  in  fact  shown  to  be  the 
case  when  comparing  the  finite  element  plastic  strain  energy  density  Wp  with  that  of  We, 
as  can  be  seen  in  Figures  6.21  and  6.22,  which  are  contours  of  the  difference  between  the 
two  energies.  As  in  the  case  of  the  plane  stress  condition,  Wp  for  the  plane  strain 
condition  is  shown  to  be  greater  than  that  of  We.  However,  both  the  amount  of  error  and 
region  of  significant  error  (greater  than  1%),  is  much  improved  over  the  plane  stress 
condition.  At  the  notch  root  itself,  the  difference  between  the  two  calculations  was  only 
7.76%  for  the  narrow  plate  and  10.6%  for  the  wide  plate. 
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Figure  6.17.  Plastic  Strain  Energy  Density  Wpfor  Narrow  Plate  in  Plane 

Strain  at  Load  Set  21. 
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Figure  6.18.  Plastic  Strain  Energy  Density  Wpfor  Wide  Plate  in  Plane  Strain 

at  Load  Set  21. 
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Figure  6.19.  Strain  Energy  Density  along  x-axisfor  Narrow  Plate  in  Plane  Strain. 


800  n 


0.0 


-wP- 


WP 


Strain  Energy  Density  for 
Elasto-Plastic  Material 

Strain  Energy  Density  for 
Elastic  Material. 


Load  Step  13j 


-[Load  Step  5 


1.0  2.0 

Distance  from  Notch  Root  (in) 


3.0 


4.0 


Figure  6.20.  Strain  Energy  Density  along  x-axisfor  Wide  Plate  in  Plane  Stress. 
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Figure  6.21.  Difference  Between  W[}  and  Wefor  Narrow  Plate  in  Plane  Strain  at 

Load  Step  21. 
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Figure  6.22.  Difference  Between  Wp  and  We  for  Wide  Plate  in  Plane  Strain  at 

Load  Step  21. 
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VII.      STRESS  AND  STRAIN  CALCULATIONS 

A.         PLANE  STRESS  CONDITION 

1.  Finite  Element  Method  Results 

The  stress  distributions  for  plane  stress  for  both  geometries  are  shown  in 
Figures  7. 1  through  7.4.  These  plots  depict  either  ay  or  ox  along  the  minimum  cross 
section  (  y  =  0  )  as  a  function  of  the  nominal  loading.  As  stated  in  Chapter  IV,  the 
maximum  Gy  value  progresses  inward  from  the  notch  root  as  yielding  increases.  This  can 
be  attributed  to  the  ax  component,  which  starts  from  zero  at  the  notch  root  and  rapidly 
approaches  its  maximum  value  inward  from  the  notch  root.  This  increase  in  the  ax 
component  results  in  a  higher  allowable  oy  than  at  the  notch  root  before  yielding  occurs. 
If  one  considers  the  limiting  case  of  a  perfectly  plastic  material;  then  at  the  notch  root, 
ay,  would  remain  at  the  value  of  the  yield  stress  after  initial  yielding.  However,  at  some 
point  inward  from  the  notch  root,  ay  would  be  higher  than  the  yield  stress  due  to  the  ax 
component.  The  other  limiting  case  is  an  elastic  only  material,  where  Gy  would  linearly 
increase  with  loading.  As  can  be  seen  from  the  stress  distributions,  the  results  of  the 
work-hardening  material  falls  somewhere  between  these  two  extremes. 
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Figure  7.1.  Distribution  of  Cx  Along  Minimum  Cross  Section  for  Narrow  Plate  in 

Plane  Stress. 
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Figure  7.2.  Distribution  of  Gv  Along  Minimum  Cross  Section  for  Narrow  Plate  in 
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Figure  7.3.  Distribution  of  Gx  Along  Minimum  Cross  Section  for  Wide  Plate  in 

Plane  Stress. 
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Figure  7.4.  Distribution  ofcyAlong  Minimum  Cross  Section  for  Wide  Plate  in 

Plane  Stress. 
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2.  Notch  Root  Stress  and  Strain  Calculations  and  Comparisons 

In  addition  to  comparing  the  plastic  strain  energy  density  Wp  to  the  predicted 
strain  energy  density  We,  comparisons  were  also  made  between  the  finite  element  stress 
and  strain  results  with  those  predicted  by  the  Glinka  and  Neuber  methods.  The  form  of 
the  Glinka  method  used  was  Equation  2.18,  while  for  the  Neuber  method,  Equation  2.4 
was  used.  If  the  nominal  stresses  are  high  enough  then  S  *  Ee,  and  Equations  2. 18  and 
2.4  are  not  valid.  They  may  be  modified,  however,  by  using  the  Ramberg-Osgood 
equation  instead  of  Hooke's  law  to  determine  the  nominal  strains  [Ref.  1,  page  138].  For 
example,  the  Neuber's  Method,  from  Equation  2.2  would  become: 


K.2  S 


S 

— h 
E 


kK) 


E         {K 


(7.1) 


Likewise,  the  Glinka  method,  from  Equation  2.12,  would  result  in: 
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(7.2) 


However,  when  these  modified  forms  were  applied  to  the  loading  levels  analyzed  in  this 
thesis,  the  amount  of  improvement  was  minimal  to  none.  For  example,  when  this 
adjustment  was  calculated  for  the  final  load  level  in  the  case  of  the  narrow  plate  under 
plane  stress  conditions,  the  error  in  notch  root  strains  at  the  final  load  step  decreased  from 
-1 1.85%  to  -1 1.66%  for  the  Glinka  method,  and  for  the  Neuber  method  the  error  actually 
increased  from  14.71%  to   14.88%. 

In  Chapter  VI,  the  strain  energy  density  calculations  showed  that  the  predicted 
strain  energy  density  based  on  elastic  material  properties  (We)  was  less  than  the  actual 
strain  energy  density  (Wp).  From  this  comparison  of  strain  energy  densities,  it  was  known 
that  the  Glinka  method  would  under  predict  the  stresses  and  strains  at  the  notch  root.  The 
stress  and  strains  were  also  predicted  based  on  the  Neuber  method.  As  was  stated  in 
Chapter  II,  the  Neuber  method  has  been  shown  to  overestimate  the  stresses  and  strains. 
This  was  also  true  for  all  of  the  configurations  analyzed  in  this  thesis.  These  results  of 
the  Glinka  method  under  predicting  the  stresses  and  strains  and  the  Neuber  method  over 
predicting  the  stresses  and  strains  are  shown  in  Figures  7.5  and  7.6. 
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Figure  7.5.  Stress/Strain  Concentration  Factors  for  Narrow  Plate  in  Plane  Stress. 
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Figure  7.6.  Stress/Strain  Concentration  Factors  for  Narrow  Plate  in  Plane  Stress. 
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These  plots  show  the  stress  and  strain  concentration  factor  as  a  function  of  the 
nominal  load  for  the  FEM  analysis,  the  Glinka  method  and  the  Neuber  method. 
Additionally,  a  stress  based  on  the  average  of  the  Glinka  and  Neuber  determined  stresses, 
and  a  strain  determined  from  the  Ramberg-Osgood  equation  and  this  average  stress  was 
determined  for  all  configurations.  It  should  be  noted  that  whenever  strain  values  are 
based  on  an  average  of  the  Glinka  and  Neuber  methods,  this  implies  that  it  is  based  on  the 
average  of  the  calculated  stresses,  and  not  the  average  of  the  strains.  The  amount  of  error 
for  the  stress  and  strain  predictions  for  both  methods  is  shown  in  Figures  7.7  through 
7. 10.  Note  that  the  strain  is  more  sensitive  than  the  stress  for  both  calculations.  This  is 
self-evident  from  the  fact  that  for  the  uniaxial  stress-strain  curve  beyond  the  yield  point, 
strain  is  highly  sensitive  to  changes  in  stress.  This  sensitivity  is  plotted  in  Figures  7.11 
and  7.12,  which  shows  the  percent  error  of  the  Glinka  method  predictions  with  respect  to 
the  percent  difference  in  strain  energy  density.  Note  that  this  sensitivity  is  also  dependent 
on  loading  condition  as  it  relates  to  the  current  stress-strain  relationship.  For  example,  as 
the  loading  increases,  the  error  in  stress  prediction  appears  to  asymptote  to  a  single  value 
after  an  initial  increase,  while  the  error  in  strain  is  almost  linear  with  respect  to  the  error 
in  strain  energy  density.  As  a  material  approaches  perfectly  plastic,  then  from  Equation 
2.6,  any  change  or  error  in  Wp  will  result  in  a  linear  change  in  e,  since  a  will  approach  a 
constant  value. 

For  the  notch  geometry  and  plane  stress  condition,  the  Glinka  and  Neuber  method 
give  an  upper  and  lower  bound  to  both  the  stress  and  strain  predictions.  Results  based  on 
the  average  of  the  stresses  of  the  Glinka  and  Neuber  method  are  in  good  agreement  to  the 
FEM  results.  It  should  be  noted  after  about  a  S/a0  ratio  of  0.5  to  0.6,  the  rate  at  which  the 
error  increases  for  the  Glinka  method  appears  to  be  constant,  while  even  though  the  error 
for  the  Neuber  method  continually  increases,  that  rate  at  which  it  increases  diminishes. 
This  results  in  the  error  for  the  average  of  these  two  methods  to  reach  a  maximum 
between  a  S/ao  ratio  of  0.54  to  0.64,  then  decrease  as  the  loading  increases.  However,  it 
does  not  appear  to  asymptote  toward  zero,  but  to  merely  change  sign  as  the  error  from  the 
Glinka  method  becomes  greater  than  that  from  the  Neuber  method. 


86 


2.0%  - 

1.5%  - 

♦    Glinka  Strain 
a    Neuber  Strain 

1.0%  - 

•    Average  of  Glinka  and  Nueber 

A 

2  0.5%  - 

u 

1— 

•s  0.0%  - 

• 

A 
A 

A     •     • 
•            *     ♦     t     ♦ 

• 

o 

u  -0.5%  - 

♦    ♦ 

-1 .0%  - 

-1.5%  - 

-2.0% 


AAAAAAAA 


A      A 


0.15 


•    • 


•    • 


♦    ♦ 


0.25 


0.35 


0.45 

S/Co 


0.55 


0.65 


0.75 


Figure  7. 7.  Error  in  Notch  Stress  <Jyfor  Glinka  and  Neuber  Method,  Narrow  Plate  in 

Plane  Stress. 
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Figure  7.8.  Error  in  Notch  Strain  £yfor  Glinka  and  Neuber  Method,  Narrow  Plate  in 

Plane  Stress. 
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Figure  7.9.  Error  in  Notch  Stress  <Jyfor  Glinka  and  Neuber  Method,  Wide  Plate  in  Plane 

Stress. 
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Figure  7.10.  Error  in  Notch  Strain  £yfor  Glinka  and  Neuber  Method,  Wide  Plate  in 

Plane  Stress. 
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Error,  Wide  Plate  in  Plane  Stress. 
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B.         PLANE  STRAIN  CONDITION 

1.  Finite  Element  Method  Results 

The  stress  distribution  for  the  plane  strain  case  are  shown  in  Figures  7.13  through 
7.18.  These  plots  depict  either  Gx,  ay  or  Gz  along  the  minimum  cross  section  (  y  =  0  )  as  a 
function  of  the  nominal  loading.  When  compared  to  the  plane  stress  condition,  the  plane 
strain  results  in  higher  axial  stresses.  This  is  due  to  the  az  stress,  which  results  in  a  higher 
hydrostatic  pressure  for  a  given  loading,  thus  reducing  the  amount  of  plastic  growth.  The 
reduced  plastic  growth  likewise  results  in  the  stresses  increasing  at  a  higher  rate  than  in 
the  plane  stress  condition.  Additionally,  the  plane  strain  condition  results  in  higher  stress 
gradients  in  the  vicinity  of  the  notch.  As  in  the  plane  stress  condition,  the  oy  stress  peaks 
at  a  point  inward  from  the  notch  root  as  plastic  growth  occurs.  For  the  plane  strain 
condition,  this  also  occurs  for  the  az  stress. 

2.  Notch  Root  Stress  and  Strain  Calculations 

For  the  plane  strain  analysis  at  the  notch  root,  Wp  is  a  function  of  ay  and  ey.  For 
the  plane  strain  condition,  however,  ey  is  a  function  of  both  ay  and  az,  and  the  strain 
energy  density  equation  does  not  reduce  to  the  simple  uniaxial  version.  To  solve  the 
plane  strain  problem,  Glinka  [Ref.  4,  5,  6]  uses  the  transformation  as  suggested  by 
Dowling,  et  al.  [Ref.  19]  that  relates  the  uniaxial  stress-strain  curve  to  a  plane  strain 
stress-strain  curve.  From  this  transformed  stress-strain  curve,  ey  can  be  found  directly 
from  ay.  To  obtain  ey,  ey  and  ay  are  related  to  the  uniaxial  stress-strain  curve  as  shown 
below: 

q,=  i  ,       ey=  /  i  (7-3) 


y 


.  2a 

where  |i  =  — 


a 


l  +  £^ 
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The  Ramgerd-Osgood  relationship  can  then  be  fitted  to  these  transformed  stresses  and 
strains,  and  new  parameters  are  determined: 


(<3     V/n" 


K 


where      E'  =  j rr  (  7.4  ) 

(1-v  j 


For  the  7075-T6  material  parameters,  the  plane  strain  transformed  parameters  are 
shown  below  in  Table  7-1.  Once  these  values  are  found,  the  stresses  at  the  notch  root  can 
then  be  found  by  substituting  all  the  normal  coefficients  with  the  transformed  coefficients 
in  Equation  2.18  and  Equation  2.4,  then  using  Equation  7.4  to  determine  the  strains. 


Coefficient  Value 


E*  (ksi) 

11,222 

K*  (ksi) 

106.77 

* 
n 

0.0541 

Table  7-1.  Plane  Strain  Transformed  Material  Coefficients. 

Since  the  results  of  the  strain  energy  density  comparison  was  much  improved  for 
the  plain  strain  condition,  it  was  expected  that  the  stress  calculations  using  the  Glinka 
method  would  also  be  improved.  This  was  indeed  the  case,  as  is  shown  in  Figures  7.19 
and  7.20,,  which  compares  the  stress  concentration  factor  as  a  function  of  the  nominal 
loading.  Note  that  again  an  average  value  of  the  two  methods  was  calculated,  and  plotted 
for  comparison.  Another  plot  of  the  strain  values  is  shown  in  Figures  7.21  and  7.22. 
Here,  K[S  is  plotted  versus  the  calculated  notch  root  strain  Ey.  Overlaid  on  this  plot  is  the 
uniaxial  stress-strain  curve,  the  transformed  plane  strain  stress-strain  curve,  and  the  notch 
stress  versus  notch  strain  results.  Not  only  does  this  show  that  the  Glinka  method  gives 
better  results  than  the  Neuber  method,  it  also  shows  that  the  FEM  notch  root  stresses  and 
strains  follow  the  transformed  stress-strain  curve.  Therefore,  one  can  conclude  that  the 
transformation  used  to  obtain  the  plane  strain  results  is  valid  at  the  notch  root.  The  error 
in  both  the  stress  and  strain  calculations  for  both  geometries  is  shown  in  Figures  7.23 
through  7.26.  These  show  that  the  Glinka  method  provides  results  three  times  more 
accurate  than  the  Neuber  method.  Additionally,  strains  calculated  based  on  the  average 
stresses  of  the  two  methods  give  slightly  more  accurate  values  than  the  Glinka  method 
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itself.  This  process  also  results  in  strains  that  are  slightly  higher  than  those  of  the  FEM 
analysis,  vice  those  of  the  Glinka  method,  which  are  slightly  lower  than  the  FEM 
analysis. 


92 


Stress,  x-component  (ksi) 

20  r 


S  (ksi) 


Distance  from  notch  root  (inches) 


Figure  7. 13.  Distribution  ofcx  Along  Minimum  Cross  Section  for  Narrow  Plate 

in  Plane  Strain. 
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Figure  7.14.  Distribution  ofo  Along  Minimum  Cross  Section  for  Narrow  Plate 
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Figure  7.15.  Distribution  of<5z  Along  Minimum  Cross  Section  for  Narrow  Plate  in 

Plane  Strain. 
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Figure  7.17.  Distribution  ofoy  Along  Minimum  Cross  Section  for  Wide  Plate 

Plane  Strain. 
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Figure  7. 18.  Distribution  of<5z  Along  Minimum  Cross  Section  for  Wide  Plate  i 

Plane  Strain. 
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Figure  7.19.  Stress/Strain  Concentration  Factors  for  Narrow  Plate  in  Plane  Strain. 
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Figure  7.20.  Stress/Strain  Concentration  Factors  for  Wide  Plate  in  Plane  Strain. 
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Figure  7.21.  Stress  and  Stress  Concentration  versus  Notch  Root  Strain  £yfor  Narrow 

Plate  in  Plane  Strain. 


o       FEM 

♦  Glinka  Method 
A       Neuber  Method 

7075-T6  Stress-Strain  Curve 

Transformed  Plane  Strain  Stress  Strain  Curve 

•  Notch  Root  Stress-Strain 


0.000        0.002        0.004        0.006        0.008        0.010        0.012        0.014        0.016 

£  (in/in) 

Figure  7.22.  Stress  and  Stress  Concentration  versus  Notch  Root  Strain  £yfor  Wide  Plate 

in  Plane  Strain. 
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Figure  7.23.  Error  in  Notch  Stress  Gyfor  Glinka  and  Neuber  Method,  Narrow  Plate  in 

Plane  Strain. 


2.0% 
1.5% 


©  1 .0%  -i 

c 

o 

I* 

S  0.5% 


co  r\  no/ 


0.0% 


w-0.5%  - 


■1.0% 


-1.5% 


AAA 


A      A 


AAA 


•      #     •      • 


'  ?  ♦ 


♦   ♦ 


♦  Glinka  method 
a  Neuber  method 

•  Average  of  Glinka  and  Neuber 


0.15 


0.25 


0.35 


0.45 

S/CTo 


0.55 


♦    ♦    ♦ 

0.65  0.75 


Figure  7.24.  Error  in  Notch  Stress  Oyfor  Glinka  and  Neuber  Method,  Wide  Plate  in 

Plane  Strain. 
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Figure  7.25.  Error  in  Notch  Strain  £yfor  Glinka  and  Neuber  Method,  Narrow  Plate  in 

Plane  Strain. 
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Figure  7.26.  Error  in  Notch  Strain  £yfor  Glinka  and  Neuber  Method,  Wide  Plate  in  Plane 

Strain. 
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VIII.     APPLICATION  OF  RESULTS  TO  FATIGUE  CALCULATIONS 

A.         STRAIN  LIFE  CALCULATIONS 

To  determine  the  consequences  of  the  Glinka  method  on  fatigue  life  calculations, 
the  results  obtained  in  Chapter  VII  for  the  plane  stress  condition  were  applied  to  a  strain 
life  analysis  as  given  in  Equation  8.1  [Ref.  1].  From  this  relationship,  the  number  of 
reversals  to  crack  initiation  (  Nf  )  is  based  on  the  strain  amplitude  (  Ae  ).  The  fatigue 
strength  coefficient  (  a'/),  fatigue  ductility  coefficient  (  £/),  fatigue  strength  exponent 
(  b  ),  and  fatigue  ductility  exponent  (  c  )  are  material  properties  and  are  listed  in  Table 


2        E 

9.1.  For  all  fatigue  calculations,  it  was  assumed  that  the  cyclic  stress-strain  curve 
remained  constant  (i.e.,  no  hysteresis  effects).  Although  this  does  not  accurately  depict 
the  cyclic  7075-T6  aluminum  properties,  it  does  provide  a  means  to  compare  the 
individual  methods  as  they  apply  to  fatigue  calculations  using  the  results  shown  in 
Chapter  VII.  Additionally,  the  fatigue  life  calculations  are  based  on  fully  reversed 
loading,  with  the  load  levels  as  shown  in  Table  6-2  used  as  the  amplitude  of  the 
alternating  load.  With  these  assumptions,  the  strain  values  obtained  previously  in 
Chapter  VII  were  simply  half  of  the  strain  amplitudes  for  the  cyclic  fatigue  calculations. 


Coefficient 

Value 

of  (ksi) 

191.0 

ef 

0.19 

b 

-0.126 

c 

-0.52 

Table  8-1.  Fatigue  Properties  of  7075 -T6  Aluminum. 
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B. 


RESULTS  OF  STRAIN  LIFE  FATIGUE  PREDICTIONS 


Figures  8.1  and  8.2  show  reversals  to  failure  (crack  initiation),  Nf,  as  predicted  by 
the  FEM,  Glinka  method,  and  Neuber  method.  These  figures  show  the  predicted  life  as  a 
function  of  the  far-field  loading.  The  fatigue  life  range  covers  from  over  6,000,000  cycles 
at  the  initial  load  for  the  narrow  plate,  which  can  be  considered  an  infinite  life,  to  less 
than  1,000  cycles  at  the  high  loads.  Figures  8.3  and  8.4  show  the  error  in  fatigue 
predictions  based  on  the  Glinka  and  Neuber  results  as  compared  to  the  FEM  analysis.  It 
is  not  surprising  that  the  greatest  error  occurs  at  the  higher  loading,  which  corresponds  to 
where  the  highest  differences  in  strain  calculations  occur.  Additionally,  the  Glinka 
method  produces  greater  errors  at  the  higher  loads  than  the  Neuber  method.  While  the 
Neuber  method  appears  to  have  a  maximum  error  just  under  40%  between  an  applied 
load  of  ±30  ksi  to±35  ksi,  the  Glinka  method  continually  gets  worse,  obtaining  a  91% 
error  for  the  wide  plate  at  fully  reversed  load  of  ±39.6  ksi. 
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Figure  8.1.  Reversals  to  Failure  for  Narrow  Plate  in  Plane  Stress. 
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Figure  8.2.  Reversals  to  Failure  for  Wide  Plate  in  Plane  Stress. 
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Figure  8.3.  Error  in  Fatigue  Life  for  Narrow  Plate  in  Plane  Stress. 
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Figure  8.4.  Error  in  Fatigue  Life  for  Wide  Plate  in  Plane  Stress. 

The  comparisons  of  the  Glinka  method  versus  the  Neuber  method  applied  to 
strain  life  predictions  show  that  although  the  error  in  strain  calculations  are  of  the  same 
magnitude,  underestimating  the  stresses  and  strains  results  in  a  greater  error  when 
calculating  the  fatigue  life.  Additionally,  not  only  does  this  result  in  a  greater  error  in 
fatigue  calculations,  but  it  errs  on  the  high  side  of  these  life  calculations.  This  could  have 
alarming  consequences  when  applied  to  safety  critical  parts.  It  should  be  noted  that  the 
strain-life  calculations  give  only  a  model  of  when  the  actual  crack  initiation  will  occur. 
These  comparisons  do,  however,  provide  an  accurate  picture  of  the  trends  in  using  either 
the  Neuber  or  Glinka  methods  in  a  popular  model  that  is  used  to  make  fatigue  life 
predictions. 
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IX.      CONCLUSIONS 

This  study  has  examined  the  proposal  by  Glinka  that  the  strain  energy  density  at 
the  notch  root  is  the  same  regardless  if  the  material  is  elastic  or  elasto-plastic.  A  detailed 
comparison  of  the  two  strain  energy  densities  was  performed  not  only  at  the  notch  root, 
but  throughout  the  field  of  symmetrical,  semi-circular  double  notched  plates.  These 
comparisons  were  made  for  both  plane  stress  and  plain  strain  conditions.  Strain  energy 
density  was  calculated  based  on  finite  element  analyses  that  had  been  rigorously  tested 
with  analytic  solutions  and  experimental  data.  The  strain  energy  density  was  numerically 
integrated,  applying  2 1  load  steps  to  reach  a  nominal  stress  equal  to  three  fourths  of  the 
yield  stress  for  each  configuration. 

In  addition  to  strain  energy  density  calculations,  stress  and  strain  calculations 
based  on  Glinka's  strain  energy  density  proposal  and  the  Neuber  method  were  performed, 
and  compared  with  the  finite  element  method  results.  To  calculate  the  notch  root  stresses 
and  strains,  a  transformation  to  a  plane  strain  stress-strain  relationship  was  performed. 
For  the  plane  stress  condition,  the  strain  results  were  applied  to  fatigue  life  predictions, 
using  the  relationship  between  number  of  reversals  and  the  strain  amplitude. 

A.         FINITE  ELEMENT  CALCULATIONS 

Prior  to  calculating  the  strain  energy  density  for  the  notched  geometries,  the  finite 
element  program  I-DEAS  Master  Series    was  thoroughly  tested  by  applying  it  to 
problems  that  offered  an  analytic  solution  and  experimental  results.  The  conclusion 
concerning  the  accuracy  of  the  finite  element  program  are  listed  below: 

•  The  finite  element  program  provided  very  accurate  results  for  the  elastic 
analysis  of  an  elliptical  hole  in  an  infinite  plate  under  uniaxial  tension.  This 
included  a  comparison  of  all  stress  components  along  the  axis  of  the  plate  and 
the  edge  of  the  ellipse. 

•  The  finite  element  analysis  of  a  hole  in  an  infinite  plate  under  elasto-plastic 
loading  again  produced  exceptional  results  when  compared  to  the  analytic 
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solution.  This  included  the  both  stress  and  strain  results  in  the  plastic  and 
elastic  region. 

•  The  finite  element  analysis  of  the  narrow  plate  with  a  central  hole  produced 
good  results  when  compared  to  the  experimental  data  of  Theocaris  and 
Marketos  [Ref.  17],  especially  at  the  hole  edge.  When  equilibrium  at  the 
minimum  section  was  tested,  the  finite  element  analysis  gave  exceptional 
results,  while  the  experimental  data  was  in  error  by  as  much  as  10%. 

The  results  show  that  the  finite  element  analysis  had  been  extensively  verified 
with  analytic  solutions  based  on  the  governing  equations  of  mechanics.  The  stresses  and 
strains  from  the  finite  element  program  were  input  into  an  integration  algorithm  that 
calculated  strain  energy  density.  Since  the  stresses  and  strains  have  been  shown  to  be 
highly  accurate,  it  can  be  inferred  that  the  strain  energy  density  calculations  are  also 
accurate. 

B.         STRAIN  ENERGY  DENSITY  CALCULATIONS 

The  strain  energy  density  based  on  an  elastic  material  and  an  elasto-plastic 
material  was  calculated  at  every  second  load  step  for  a  total  of  twenty  one  loads.  A 
detailed  comparisons  of  the  two  strain  energy  densities  reveal  that: 

•  The  strain  energy  density  in  the  vicinity  of  the  notch  root  based  on  elasto- 
plastic  material  properties  is  higher  than  the  strain  energy  density  assuming 
elastic  only  properties. 

•  The  plane  strain  condition  results  in  better  agreement  between  the  two  strain 
energy  densities  than  the  plane  stress  condition. 

•  At  the  higher  loads,  the  greatest  deviation  between  the  two  energies  occurs 
slightly  inward  from  the  notch  root,  vice  at  the  notch  root  itself. 

•  The  distribution  of  the  two  strain  energy  densities  along  the  minimum  cross 
section  not  only  differ  in  magnitude,  but  also  in  shape. 
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C.  NOTCH  ROOT  STRESS  AND  STRAIN  CALCULATIONS 

The  notch  root  stress  and  strain  values  were  calculated  based  on  the  Glinka 
proposal  at  the  notch  root,  and  compared  to  the  finite  element  data  and  the  Neuber 
method  of  estimating  notch  root  stresses  and  strains.  The  comparison  of  these  two 
methods  revealed: 

•  The  Glinka  strain  energy  density  method  under  estimates  the  stresses  and 
strains,  while  the  Neuber  method  overestimates  the  stresses  and  strains. 

•  For  the  plane  stress  condition,  the  two  methods  appeared  to  give  an  upper  and 
lower  bound.  Taking  the  average  of  the  stresses  from  each  method  and 
determining  the  strains  from  this  value  gave  very  good  results. 

•  The  Glinka  method,  while  under  estimating  the  stress  and  strain  values, 
produces  results  two  to  three  times  better  than  the  Neuber  method  for  plane 
strain  conditions. 

D.  IMPACT  ON  FATIGUE  LIFE  PREDICTIONS 

Fatigue  life  calculations  were  made  based  on  the  Glinka  method,  the  Neuber 
method,  and  the  average  of  the  stress  values  of  these  two  methods,  and  compared  to 
fatigue  life  calculations  based  on  the  FEM  results  for  the  plane  stress  condition.  The 
impact  of  the  two  methods  is  summarized  below: 

•  Since  the  Glinka  method  under  estimated  the  stresses  and  strains,  it 
overestimated  the  predicted  life.  However,  the  amount  of  error  steadily 
increased  with  the  applied  load,  growing  up  to  90%  for  the  wide  plate 
geometry  at  its  final  load  level.  For  the  Neuber  method,  the  error  appears  to 
reach  a  maximum  of  40%  for  a  lower  cyclic  load  and  improves  slightly  as  the 
load  increases. 

•  Estimating  the  fatigue  life  by  using  an  average  of  the  stresses  of  the  two 
methods  produce  results  that  were  ±10%  from  those  based  on  the  finite 
element  results. 
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Even  though  the  fatigue  calculations  were  based  on  finite  element  data,  with  no 
comparisons  made  to  experimental  data,  a  valid  comparison  between  the  two  methods 
was  made.  Since  the  Glinka  method  under  estimates  the  stresses  and  strains,  the  fatigue 
life  based  on  this  method  will  be  greater  than  the  actual  fatigue  life.  This  would  result  in 
parts  failing  prior  to  their  expected  life  cycle,  and  would  be  very  detrimental  unless  safety 
factors  were  built  into  the  structure  and  the  calculations. 

E.         RECOMMENDATIONS 

This  study  thoroughly  analyzed  the  Glinka  strain  energy  density  proposal  for 
stress  concentration  calculations  at  the  notch  root,  and  compared  the  findings  with  results 
from  the  Neuber  method.  For  the  geometries  studied,  it  was  shown  that  the  Glinka 
method  under-predicts  the  stresses  and  strains,  while  the  Neuber  method  over-predicts  the 
stresses  and  strains.  Further  comparisons  should  be  made  not  only  to  different 
geometries,  but  also  with  different  material  properties.  A  study  should  also  be  performed 
by  comparing  actual  fatigue  data  with  predictions  based  on  the  two  methods.  Taking  an 
average  value  of  the  two  methods  appears  to  give  good  results,  and  as  in  the  case  of  either 
the  Glinka  or  Neuber  method,  can  be  a  means  of  quickly  and  efficiently  computing  notch 
root  stresses  and  strains.  This  mean  value  of  the  two  methods  should  also  be  compared 
with  different  geometries  and  material  properties,  along  with  comparisons  to  actual 
fatigue  data. 
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